MODULE SPECDATA

USE PRECISION_PARAMETERS
IMPLICIT NONE (TYPE,EXTERNAL)
INTEGER ::  I,J
INTEGER, PARAMETER :: NWATERK=183
REAL(EB) :: CPLXREF_WATER(NWATERK,2)
INTEGER, PARAMETER :: NFUELK=94
REAL(EB) :: CPLXREF_FUEL(NFUELK,3)

! HALE, G.M. AND QUERRY, M.R., "OPTICAL CONSTANTS OF WATER IN THE 200-NM TO 200-\MU M WAVELENGTH REGION",
! APPLIED OPTICS, 12(3), 555-63 (1973)

DATA ((CPLXREF_WATER(I,J), J = 1,2), I = 1,10)  &
/ 1.0000000E-6_EB,  2.8647890E-6_EB,&
1.0200000E-6_EB,  2.1915636E-6_EB,&
1.0400000E-6_EB,  1.3241691E-6_EB,&
1.0600000E-6_EB,  1.0122254E-6_EB,&
1.0800000E-6_EB,  1.1172677E-6_EB,&
1.1000000E-6_EB,  1.4880987E-6_EB,&
1.1200000E-6_EB,  4.6345919E-6_EB,&
1.1400000E-6_EB,  5.9874090E-6_EB,&
1.1600000E-6_EB,  8.2155782E-6_EB,&
1.1800000E-6_EB,  9.7657473E-6_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =11,20)  &
/ 1.2000000E-6_EB,  9.9312684E-6_EB,&
1.2200000E-6_EB,  9.2230290E-6_EB,&
1.2400000E-6_EB,  8.6834937E-6_EB,&
1.2600000E-6_EB,  8.9238177E-6_EB,&
1.2800000E-6_EB,  9.9821980E-6_EB,&
1.3000000E-6_EB,  1.1483029E-5_EB,&
1.3200000E-6_EB,  1.1953809E-5_EB,&
1.3400000E-6_EB,  1.2614780E-5_EB,&
1.3600000E-6_EB,  2.9978533E-5_EB,&
1.3800000E-6_EB,  6.6988316E-5_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =21,30)  &
/ 1.4000000E-6_EB,  1.3803508E-4_EB,&
1.4200000E-6_EB,  2.4995602E-4_EB,&
1.4400000E-6_EB,  3.3002369E-4_EB,&
1.4600000E-6_EB,  3.2996003E-4_EB,&
1.4800000E-6_EB,  2.5003560E-4_EB,&
1.5000000E-6_EB,  2.0996516E-4_EB,&
1.5200000E-6_EB,  1.6994565E-4_EB,&
1.5400000E-6_EB,  1.4497583E-4_EB,&
1.5600000E-6_EB,  1.2004433E-4_EB,&
1.5800000E-6_EB,  9.9957388E-5_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =31,40)  &
/ 1.6000000E-6_EB,  8.5561825E-5_EB,&
1.6200000E-6_EB,  7.5028952E-5_EB,&
1.6400000E-6_EB,  6.4992643E-5_EB,&
1.6600000E-6_EB,  5.9972898E-5_EB,&
1.6800000E-6_EB,  6.0027012E-5_EB,&
1.7000000E-6_EB,  6.0065211E-5_EB,&
1.7200000E-6_EB,  6.9942231E-5_EB,&
1.7400000E-6_EB,  8.5017388E-5_EB,&
1.7600000E-6_EB,  1.0000023E-4_EB,&
1.7800000E-6_EB,  1.1501809E-4_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =41,50)  &
/ 1.8000000E-6_EB,  1.1502128E-4_EB,&
1.8200000E-6_EB,  1.3005838E-4_EB,&
1.8400000E-6_EB,  1.4993669E-4_EB,&
1.8600000E-6_EB,  2.1003200E-4_EB,&
1.8800000E-6_EB,  4.6497435E-4_EB,&
1.9000000E-6_EB,  1.0000183E-3_EB,&
1.9200000E-6_EB,  1.7500423E-3_EB,&
1.9400000E-6_EB,  1.8499391E-3_EB,&
1.9600000E-6_EB,  1.6500261E-3_EB,&
1.9800000E-6_EB,  1.4500559E-3_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =51,60)  &
/ 2.0000000E-6_EB,  1.1000790E-3_EB,&
2.0200000E-6_EB,  9.0001961E-4_EB,&
2.0400000E-6_EB,  7.3003417E-4_EB,&
2.0600000E-6_EB,  6.3998112E-4_EB,&
2.0800000E-6_EB,  5.2006742E-4_EB,&
2.1000000E-6_EB,  4.5003447E-4_EB,&
2.1200000E-6_EB,  4.0505888E-4_EB,&
2.1400000E-6_EB,  3.4995785E-4_EB,&
2.1600000E-6_EB,  3.2005422E-4_EB,&
2.1800000E-6_EB,  2.9994500E-4_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =61,70)  &
/ 2.2000000E-6_EB,  2.8904129E-4_EB,&
2.2200000E-6_EB,  2.8495578E-4_EB,&
2.2400000E-6_EB,  2.9500960E-4_EB,&
2.2600000E-6_EB,  3.1005293E-4_EB,&
2.2800000E-6_EB,  3.5997028E-4_EB,&
2.3000000E-6_EB,  4.0998313E-4_EB,&
2.3200000E-6_EB,  4.9496551E-4_EB,&
2.3400000E-6_EB,  5.9494505E-4_EB,&
2.3600000E-6_EB,  6.9994116E-4_EB,&
2.3800000E-6_EB,  8.2007768E-4_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =71,80)  &
/ 2.4000000E-6_EB,  9.5607557E-4_EB,&
2.4200000E-6_EB,  1.1500727E-3_EB,&
2.4400000E-6_EB,  1.2999617E-3_EB,&
2.4600000E-6_EB,  1.4999176E-3_EB,&
2.4800000E-6_EB,  1.6999912E-3_EB,&
2.5000000E-6_EB,  1.8000424E-3_EB,&
2.5200000E-6_EB,  2.0500716E-3_EB,&
2.5400000E-6_EB,  2.1999478E-3_EB,&
2.5600000E-6_EB,  2.3500946E-3_EB,&
2.5800000E-6_EB,  2.7000302E-3_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =81,90)  &
/ 2.6000000E-6_EB,  3.1699367E-3_EB,&
2.6500000E-6_EB,  6.7000889E-3_EB,&
2.7000000E-6_EB,  1.8999997E-2_EB,&
2.7500000E-6_EB,  5.9000926E-2_EB,&
2.8000000E-6_EB,  1.1500027E-1_EB,&
2.8500000E-6_EB,  1.8499960E-1_EB,&
2.9000000E-6_EB,  2.6769861E-1_EB,&
2.9500000E-6_EB,  2.9813700E-1_EB,&
3.0000000E-6_EB,  2.7215495E-1_EB,&
3.0500000E-6_EB,  2.4000020E-1_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =91,100)  &
/ 3.1000000E-6_EB,  1.9200142E-1_EB,&
3.1500000E-6_EB,  1.3500032E-1_EB,&
3.2000000E-6_EB,  9.2401540E-2_EB,&
3.2500000E-6_EB,  6.0999713E-2_EB,&
3.3000000E-6_EB,  3.6798931E-2_EB,&
3.3500000E-6_EB,  2.6099958E-2_EB,&
3.4000000E-6_EB,  1.9500046E-2_EB,&
3.4500000E-6_EB,  1.3199993E-2_EB,&
3.5000000E-6_EB,  9.4000888E-3_EB,&
3.6000000E-6_EB,  5.1500311E-3_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =101,110)  &
/ 3.7000000E-6_EB,  3.6000769E-3_EB,&
3.8000000E-6_EB,  3.4001225E-3_EB,&
3.9000000E-6_EB,  3.7999516E-3_EB,&
4.0000000E-6_EB,  4.5998962E-3_EB,&
4.1000000E-6_EB,  5.6118033E-3_EB,&
4.2000000E-6_EB,  6.8850428E-3_EB,&
4.3000000E-6_EB,  8.4519233E-3_EB,&
4.4000000E-6_EB,  1.0294142E-2_EB,&
4.5000000E-6_EB,  1.3392888E-2_EB,&
4.6000000E-6_EB,  1.4715466E-2_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =111,120)  &
/ 4.7000000E-6_EB,  1.5708593E-2_EB,&
4.8000000E-6_EB,  1.5011494E-2_EB,&
4.9000000E-6_EB,  1.3686529E-2_EB,&
5.0000000E-6_EB,  1.2414086E-2_EB,&
5.1000000E-6_EB,  1.1120156E-2_EB,&
5.2000000E-6_EB,  1.0096790E-2_EB,&
5.3000000E-6_EB,  9.7848459E-3_EB,&
5.4000000E-6_EB,  1.0313240E-2_EB,&
5.5000000E-6_EB,  1.1598416E-2_EB,&
5.6000000E-6_EB,  1.4215720E-2_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =121,130)  &
/ 5.7000000E-6_EB,  2.0320903E-2_EB,&
5.8000000E-6_EB,  3.3000777E-2_EB,&
5.9000000E-6_EB,  6.2209688E-2_EB,&
6.0000000E-6_EB,  1.0699987E-1_EB,&
6.1000000E-6_EB,  1.3101555E-1_EB,&
6.2000000E-6_EB,  8.8019050E-2_EB,&
6.3000000E-6_EB,  5.7002139E-2_EB,&
6.4000000E-6_EB,  4.4868962E-2_EB,&
6.5000000E-6_EB,  3.9207820E-2_EB,&
6.6000000E-6_EB,  3.5609327E-2_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =131,140)  &
/ 6.7000000E-6_EB,  3.3696285E-2_EB,&
6.8000000E-6_EB,  3.2684059E-2_EB,&
6.9000000E-6_EB,  3.2176355E-2_EB,&
7.0000000E-6_EB,  3.1974228E-2_EB,&
7.1000000E-6_EB,  3.1979003E-2_EB,&
7.2000000E-6_EB,  3.2085637E-2_EB,&
7.3000000E-6_EB,  3.2182721E-2_EB,&
7.4000000E-6_EB,  3.2388031E-2_EB,&
7.5000000E-6_EB,  3.2586975E-2_EB,&
7.6000000E-6_EB,  3.2779552E-2_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =141,150)  &
/ 7.7000000E-6_EB,  3.3088313E-2_EB,&
7.8000000E-6_EB,  3.3518031E-2_EB,&
7.9000000E-6_EB,  3.3884883E-2_EB,&
8.0000000E-6_EB,  3.4313806E-2_EB,&
8.2000000E-6_EB,  3.5106397E-2_EB,&
8.4000000E-6_EB,  3.6096341E-2_EB,&
8.6000000E-6_EB,  3.1001791E-2_EB,&
8.8000000E-6_EB,  3.8515496E-2_EB,&
9.0000000E-6_EB,  3.9892186E-2_EB,&
9.2000000E-6_EB,  4.1510792E-2_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =151,160)  &
/ 9.4000000E-6_EB,  4.3310835E-2_EB,&
9.6000000E-6_EB,  4.5378257E-2_EB,&
9.8000000E-6_EB,  4.7883356E-2_EB,&
1.0000000E-5_EB,  5.0770427E-2_EB,&
1.0500000E-5_EB,  6.6176625E-2_EB,&
1.1000000E-5_EB,  9.6813952E-2_EB,&
1.1500000E-5_EB,  1.4202987E-1_EB,&
1.2000000E-5_EB,  1.9900734E-1_EB,&
1.2500000E-5_EB,  2.5902467E-1_EB,&
1.3000000E-5_EB,  3.0497270E-1_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =161,170)  &
/ 1.3500000E-5_EB,  3.4302267E-1_EB,&
1.4000000E-5_EB,  3.6998750E-1_EB,&
1.4500000E-5_EB,  3.8804760E-1_EB,&
1.5000000E-5_EB,  4.0202539E-1_EB,&
1.5500000E-5_EB,  4.1394609E-1_EB,&
1.6000000E-5_EB,  4.2195159E-1_EB,&
1.6500000E-5_EB,  4.2804722E-1_EB,&
1.7000000E-5_EB,  4.2897828E-1_EB,&
1.7500000E-5_EB,  4.2906183E-1_EB,&
1.8000000E-5_EB,  4.2599412E-1_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =171,180)  &
/ 1.8500000E-5_EB,  4.2104440E-1_EB,&
1.9000000E-5_EB,  4.1397792E-1_EB,&
1.9500000E-5_EB,  4.0407849E-1_EB,&
2.0000000E-5_EB,  3.9295355E-1_EB,&
3.0000000E-5_EB,  3.2801834E-1_EB,&
3.8000000E-5_EB,  3.6105890E-1_EB,&
5.0000000E-5_EB,  5.1407047E-1_EB,&
6.0000000E-5_EB,  5.8680428E-1_EB,&
7.0000000E-5_EB,  5.7598174E-1_EB,&
8.0000000E-5_EB,  5.4685638E-1_EB /
DATA ((CPLXREF_WATER(I,J), J = 1,2), I =181,NWATERK)  &
/ 9.0000000E-5_EB,  5.3571554E-1_EB,&
1.0000000E-4_EB,  5.3237328E-1_EB,&
2.0000000E-4_EB,  5.0452117E-1_EB /

! HEPTANE PROPERTIES FROM
! L.A. DOMBROVSKY, S.S. SAZHIN, S.V. MIKHALOVSKY, R. WOOD, M.R. HEIKAL
! SPECTRAL PROPERTIES OF DIESEL FUEL PARTICLES
! FUEL, VOL. 82, NO. 1 (2003) PP. 15-22
DATA (CPLXREF_FUEL( 1,J), J=1,3) / 0.7_EB,   1.47_EB, 1.55E-07_EB /
DATA (CPLXREF_FUEL( 2,J), J=1,3) / 0.8_EB,   1.47_EB, 4.33E-07_EB /
DATA (CPLXREF_FUEL( 3,J), J=1,3) / 0.9_EB,   1.47_EB, 9.39E-07_EB /
DATA (CPLXREF_FUEL( 4,J), J=1,3) / 1.0_EB,   1.47_EB, 8.81E-07_EB /
DATA (CPLXREF_FUEL( 5,J), J=1,3) / 1.1_EB,   1.47_EB, 1.19E-06_EB /
DATA (CPLXREF_FUEL( 6,J), J=1,3) / 1.2_EB,   1.47_EB, 2.63E-06_EB /
DATA (CPLXREF_FUEL( 7,J), J=1,3) / 1.3_EB,   1.47_EB, 1.08E-05_EB /
DATA (CPLXREF_FUEL( 8,J), J=1,3) / 1.4_EB,   1.47_EB, 2.24E-05_EB /
DATA (CPLXREF_FUEL( 9,J), J=1,3) / 1.5_EB,   1.46_EB, 3.41E-05_EB /
DATA (CPLXREF_FUEL( 10,J), J=1,3) / 1.6_EB,  1.46_EB, 4.57E-05_EB /
DATA (CPLXREF_FUEL( 11,J), J=1,3) / 1.7_EB,  1.46_EB, 5.73E-05_EB /
DATA (CPLXREF_FUEL( 12,J), J=1,3) / 1.8_EB,  1.46_EB, 6.90E-05_EB /
DATA (CPLXREF_FUEL( 13,J), J=1,3) / 1.9_EB,  1.46_EB, 8.06E-05_EB /
DATA (CPLXREF_FUEL( 14,J), J=1,3) / 2.0_EB,  1.46_EB, 9.22E-05_EB /
DATA (CPLXREF_FUEL( 15,J), J=1,3) / 2.1_EB,  1.46_EB, 1.04E-04_EB /
DATA (CPLXREF_FUEL( 16,J), J=1,3) / 2.2_EB,  1.46_EB, 1.15E-04_EB /
DATA (CPLXREF_FUEL( 17,J), J=1,3) / 2.3_EB,  1.46_EB, 1.27E-04_EB /
DATA (CPLXREF_FUEL( 18,J), J=1,3) / 2.4_EB,  1.45_EB, 1.39E-04_EB /
DATA (CPLXREF_FUEL( 19,J), J=1,3) / 2.5_EB,  1.45_EB, 1.50E-04_EB /
DATA (CPLXREF_FUEL( 20,J), J=1,3) / 2.6_EB,  1.45_EB, 1.62E-04_EB /
DATA (CPLXREF_FUEL( 21,J), J=1,3) / 2.7_EB,  1.45_EB, 1.11E-04_EB /
DATA (CPLXREF_FUEL( 22,J), J=1,3) / 2.8_EB,  1.44_EB, 5.92E-05_EB /
DATA (CPLXREF_FUEL( 23,J), J=1,3) / 2.9_EB,  1.44_EB, 7.45E-05_EB /
DATA (CPLXREF_FUEL( 24,J), J=1,3) / 3.0_EB,  1.43_EB, 9.72E-05_EB /
DATA (CPLXREF_FUEL( 25,J), J=1,3) / 3.1_EB,  1.42_EB, 3.12E-04_EB /
DATA (CPLXREF_FUEL( 26,J), J=1,3) / 3.2_EB,  1.40_EB, 6.09E-04_EB /
DATA (CPLXREF_FUEL( 27,J), J=1,3) / 3.3_EB,  1.17_EB, 5.72E-02_EB /
DATA (CPLXREF_FUEL( 28,J), J=1,3) / 3.4_EB,  1.39_EB, 1.20E-01_EB /
DATA (CPLXREF_FUEL( 29,J), J=1,3) / 3.5_EB,  1.45_EB, 8.24E-02_EB /
DATA (CPLXREF_FUEL( 30,J), J=1,3) / 3.6_EB,  1.51_EB, 1.63E-03_EB /
DATA (CPLXREF_FUEL( 31,J), J=1,3) / 3.7_EB,  1.64_EB, 1.33E-03_EB /
DATA (CPLXREF_FUEL( 32,J), J=1,3) / 3.8_EB,  1.56_EB, 1.02E-03_EB /
DATA (CPLXREF_FUEL( 33,J), J=1,3) / 3.9_EB,  1.52_EB, 6.36E-04_EB /
DATA (CPLXREF_FUEL( 34,J), J=1,3) / 4.0_EB,  1.50_EB, 2.51E-04_EB /
DATA (CPLXREF_FUEL( 35,J), J=1,3) / 4.1_EB,  1.50_EB, 2.59E-04_EB /
DATA (CPLXREF_FUEL( 36,J), J=1,3) / 4.2_EB,  1.49_EB, 3.10E-04_EB /
DATA (CPLXREF_FUEL( 37,J), J=1,3) / 4.3_EB,  1.49_EB, 2.60E-04_EB /
DATA (CPLXREF_FUEL( 38,J), J=1,3) / 4.4_EB,  1.48_EB, 2.11E-04_EB /
DATA (CPLXREF_FUEL( 39,J), J=1,3) / 4.5_EB,  1.48_EB, 1.98E-04_EB /
DATA (CPLXREF_FUEL( 40,J), J=1,3) / 4.6_EB,  1.47_EB, 1.90E-04_EB /
DATA (CPLXREF_FUEL( 41,J), J=1,3) / 4.7_EB,  1.47_EB, 1.58E-04_EB /
DATA (CPLXREF_FUEL( 42,J), J=1,3) / 4.8_EB,  1.47_EB, 1.21E-04_EB /
DATA (CPLXREF_FUEL( 43,J), J=1,3) / 4.9_EB,  1.47_EB, 8.32E-05_EB /
DATA (CPLXREF_FUEL( 44,J), J=1,3) / 5.0_EB,  1.46_EB, 4.58E-05_EB /
DATA (CPLXREF_FUEL( 45,J), J=1,3) / 5.1_EB,  1.46_EB, 6.37E-05_EB /
DATA (CPLXREF_FUEL( 46,J), J=1,3) / 5.2_EB,  1.46_EB, 7.78E-05_EB /
DATA (CPLXREF_FUEL( 47,J), J=1,3) / 5.3_EB,  1.46_EB, 7.66E-05_EB /
DATA (CPLXREF_FUEL( 48,J), J=1,3) / 5.4_EB,  1.45_EB, 7.55E-05_EB /
DATA (CPLXREF_FUEL( 49,J), J=1,3) / 5.5_EB,  1.45_EB, 9.62E-05_EB /
DATA (CPLXREF_FUEL( 50,J), J=1,3) / 5.6_EB,  1.45_EB, 1.17E-04_EB /
DATA (CPLXREF_FUEL( 51,J), J=1,3) / 5.7_EB,  1.45_EB, 1.63E-04_EB /
DATA (CPLXREF_FUEL( 52,J), J=1,3) / 5.8_EB,  1.45_EB, 2.16E-04_EB /
DATA (CPLXREF_FUEL( 53,J), J=1,3) / 5.9_EB,  1.44_EB, 2.68E-04_EB /
DATA (CPLXREF_FUEL( 54,J), J=1,3) / 6.0_EB,  1.44_EB, 3.21E-04_EB /
DATA (CPLXREF_FUEL( 55,J), J=1,3) / 6.1_EB,  1.44_EB, 3.74E-04_EB /
DATA (CPLXREF_FUEL( 56,J), J=1,3) / 6.2_EB,  1.44_EB, 4.36E-04_EB /
DATA (CPLXREF_FUEL( 57,J), J=1,3) / 6.3_EB,  1.43_EB, 5.87E-04_EB /
DATA (CPLXREF_FUEL( 58,J), J=1,3) / 6.4_EB,  1.43_EB, 7.38E-04_EB /
DATA (CPLXREF_FUEL( 59,J), J=1,3) / 6.5_EB,  1.42_EB, 1.35E-03_EB /
DATA (CPLXREF_FUEL( 60,J), J=1,3) / 6.6_EB,  1.41_EB, 6.12E-03_EB /
DATA (CPLXREF_FUEL( 61,J), J=1,3) / 6.7_EB,  1.39_EB, 2.06E-02_EB /
DATA (CPLXREF_FUEL( 62,J), J=1,3) / 6.8_EB,  1.35_EB, 3.51E-02_EB /
DATA (CPLXREF_FUEL( 63,J), J=1,3) / 6.9_EB,  1.37_EB, 2.29E-02_EB /
DATA (CPLXREF_FUEL( 64,J), J=1,3) / 7.0_EB,  1.48_EB, 3.99E-03_EB /
DATA (CPLXREF_FUEL( 65,J), J=1,3) / 7.1_EB,  1.53_EB, 3.24E-03_EB /
DATA (CPLXREF_FUEL( 66,J), J=1,3) / 7.2_EB,  1.48_EB, 2.61E-03_EB /
DATA (CPLXREF_FUEL( 67,J), J=1,3) / 7.3_EB,  1.43_EB, 2.97E-03_EB /
DATA (CPLXREF_FUEL( 68,J), J=1,3) / 7.4_EB,  1.47_EB, 3.33E-03_EB /
DATA (CPLXREF_FUEL( 69,J), J=1,3) / 7.5_EB,  1.50_EB, 2.80E-03_EB /
DATA (CPLXREF_FUEL( 70,J), J=1,3) / 7.6_EB,  1.49_EB, 2.14E-03_EB /
DATA (CPLXREF_FUEL( 71,J), J=1,3) / 7.7_EB,  1.48_EB, 2.28E-03_EB /
DATA (CPLXREF_FUEL( 72,J), J=1,3) / 7.8_EB,  1.47_EB, 2.41E-03_EB /
DATA (CPLXREF_FUEL( 73,J), J=1,3) / 7.9_EB,  1.47_EB, 1.63E-03_EB /
DATA (CPLXREF_FUEL( 74,J), J=1,3) / 8.0_EB,  1.47_EB, 9.84E-04_EB /
DATA (CPLXREF_FUEL( 75,J), J=1,3) / 8.1_EB,  1.46_EB, 9.03E-04_EB /
DATA (CPLXREF_FUEL( 76,J), J=1,3) / 8.2_EB,  1.46_EB, 8.19E-04_EB /
DATA (CPLXREF_FUEL( 77,J), J=1,3) / 8.3_EB,  1.46_EB, 7.13E-04_EB /
DATA (CPLXREF_FUEL( 78,J), J=1,3) / 8.4_EB,  1.46_EB, 6.07E-04_EB /
DATA (CPLXREF_FUEL( 79,J), J=1,3) / 8.5_EB,  1.45_EB, 8.54E-04_EB /
DATA (CPLXREF_FUEL( 80,J), J=1,3) / 8.6_EB,  1.45_EB, 1.10E-03_EB /
DATA (CPLXREF_FUEL( 81,J), J=1,3) / 8.7_EB,  1.45_EB, 1.76E-03_EB /
DATA (CPLXREF_FUEL( 82,J), J=1,3) / 8.8_EB,  1.45_EB, 2.41E-03_EB /
DATA (CPLXREF_FUEL( 83,J), J=1,3) / 8.9_EB,  1.45_EB, 1.58E-03_EB /
DATA (CPLXREF_FUEL( 84,J), J=1,3) / 9.0_EB,  1.45_EB, 5.71E-04_EB /
DATA (CPLXREF_FUEL( 85,J), J=1,3) / 9.1_EB,  1.45_EB, 8.53E-04_EB /
DATA (CPLXREF_FUEL( 86,J), J=1,3) / 9.2_EB,  1.45_EB, 1.28E-03_EB /
DATA (CPLXREF_FUEL( 87,J), J=1,3) / 9.3_EB,  1.45_EB, 1.63E-03_EB /
DATA (CPLXREF_FUEL( 88,J), J=1,3) / 9.4_EB,  1.45_EB, 1.98E-03_EB /
DATA (CPLXREF_FUEL( 89,J), J=1,3) / 9.5_EB,  1.45_EB, 1.69E-03_EB /
DATA (CPLXREF_FUEL( 90,J), J=1,3) / 9.6_EB,  1.45_EB, 1.24E-03_EB /
DATA (CPLXREF_FUEL( 91,J), J=1,3) / 9.7_EB,  1.45_EB, 1.43E-03_EB /
DATA (CPLXREF_FUEL( 92,J), J=1,3) / 9.8_EB,  1.45_EB, 1.69E-03_EB /
DATA (CPLXREF_FUEL( 93,J), J=1,3) / 9.9_EB,  1.45_EB, 1.26E-03_EB /
DATA (CPLXREF_FUEL( 94,J), J=1,3) / 10.0_EB, 1.45_EB, 8.24E-04_EB /

END MODULE SPECDATA

MODULE WSGG_ARRAYS
USE PRECISION_PARAMETERS

REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: WSGG_B1_ARRAY,WSGG_B2_ARRAY,WSGG_D_ARRAY
REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: WSGG_C_ARRAY
REAL(EB), ALLOCATABLE, DIMENSION(:) :: WSGG_KAPPAP1_ARRAY,WSGG_KAPPAP2_ARRAY

END MODULE WSGG_ARRAYS

MODULE MIEV

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS, ONLY: NUMBER_SPECTRAL_BANDS,NUMBER_RADIATION_ANGLES,LU_ERR
USE SPECDATA, ONLY: CPLXREF_WATER, NWATERK, CPLXREF_FUEL, NFUELK
USE RADCAL_CALC, ONLY: PLANCK
USE RADCONS
USE MEMORY_FUNCTIONS, ONLY : CHKMEMERR
USE MATH_FUNCTIONS, ONLY : INTERPOLATE1D

IMPLICIT NONE (TYPE,EXTERNAL)
REAL(EB), ALLOCATABLE :: RDMIE(:),  LMBDMIE(:),LMBDWGHT(:),REAL_REF_INDX(:),CMPLX_REF_INDX(:)
REAL(EB), ALLOCATABLE :: QSCA(:,:), QABS(:,:),  CHI_F(:,:)

PRIVATE
PUBLIC MEAN_CROSS_SECTIONS

CONTAINS


!> \brief Calculates the mean scattering and absorption coefficients for each radiation band and particle size group
!> \param PARTICLE_CLASS Index of the Lagrangian particle class
!> \param SMIX_INDEX Index of the species mixture

SUBROUTINE MEAN_CROSS_SECTIONS(PARTICLE_CLASS,SMIX_INDEX)

USE TYPES, ONLY: LAGRANGIAN_PARTICLE_CLASS_TYPE, LAGRANGIAN_PARTICLE_CLASS, TABLES_TYPE, TABLES, SPECIES_MIXTURE, SPECIES, &
                 SPECIES_MIXTURE_TYPE
USE GLOBAL_CONSTANTS, ONLY: H2O_INDEX
USE MATH_FUNCTIONS, ONLY : INTERPOLATE1D

INTEGER, INTENT(IN), OPTIONAL :: PARTICLE_CLASS,SMIX_INDEX
INTEGER  :: NSB,I,J,IBND,IZERO,NX,NLAMBDALOW(1),NLAMBDAHIGH(1),ND,NL_DATA
REAL(EB) :: RMMAX,RMMIN,RDTMP,IB,IBSUM,AVAL,BVAL,ASUM,BSUM,B_WIEN
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: WQABS,WQSCA
REAL(EB), ALLOCATABLE, DIMENSION(:) :: R50
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE (SPECIES_MIXTURE_TYPE), POINTER :: SM2
TYPE (TABLES_TYPE),  POINTER :: TA

IF (PRESENT(PARTICLE_CLASS)) THEN
   LPC => LAGRANGIAN_PARTICLE_CLASS(PARTICLE_CLASS)
ELSE
   SM2 => SPECIES_MIXTURE(SMIX_INDEX)
ENDIF
NSB = NUMBER_SPECTRAL_BANDS

! PHYSICAL PARAMETERS
! MINIMUM MEAN RADIUS (M)
RMMIN = 0.5_EB*1.E-6*MIE_MINIMUM_DIAMETER
IF (RMMIN < TWO_EPSILON_EB) RMMIN = 0.5E-6_EB

! MAXIMUM MEAN RADIUS (M)
RMMAX = 0.5_EB*1.E-6*MIE_MAXIMUM_DIAMETER
IF (RMMAX < TWO_EPSILON_EB) THEN
   IF (PRESENT(PARTICLE_CLASS)) THEN
      RMMAX = 0.5_EB*LPC%DIAMETER
   ELSE
      RMMAX = 0.5_EB*SM2%MEAN_DIAMETER
   ENDIF
   ! ALLOW INCREASE OF THE MEAN RADIUS
   RMMAX = 1.5_EB*RMMAX
ENDIF

! OTHER CONSTANTS

B_WIEN = 2.8977685E-3_EB

! CALCULATE PARAMETERS OF THE PARTICLE GROUP LOOKUP TABLE

DGROUP_A = (LOG(RMMAX)-LOG(RMMIN))/(MIE_NDG-1)
DGROUP_B = LOG(RMMAX)-DGROUP_A*MIE_NDG

! GENERATE THE PARTICLE RADII FOR MIE TABLE (MICRONS)

RDTMP = 0.5_EB*RMMIN
NX = 0
DO WHILE (RDTMP < 2._EB*RMMAX)
   NX = NX + 1
   RDTMP = RDTMP + MIN(1.0E6_EB,0.2_EB*RDTMP**(1._EB))
ENDDO
NRDMIE = NX

ALLOCATE(RDMIE(1:NRDMIE),STAT=IZERO)
CALL CHKMEMERR('MIEV','RDMIE',IZERO)

RDTMP = 0.5_EB*RMMIN
RDMIE(1) = RDTMP
DO NX = 2, NRDMIE
   RDTMP = RDTMP + MIN(3.0E6_EB,0.2_EB*RDTMP**(1.0_EB))
   RDMIE(NX) = RDTMP
ENDDO

! RADIATIVE PROPERTIES

NLMBDMIE = NWATERK
IF (PRESENT(PARTICLE_CLASS)) THEN
   IF (LPC%RADIATIVE_PROPERTY_INDEX>0) THEN
      TA => TABLES(LPC%RADIATIVE_PROPERTY_INDEX)
      NL_DATA = TA%NUMBER_ROWS
   ELSEIF (LPC%Y_INDEX==H2O_INDEX) THEN
      NL_DATA = NLMBDMIE
   ELSEIF (LPC%FUEL) THEN
      NL_DATA = NFUELK
   ELSE
      NL_DATA = 1
   ENDIF
ELSE
   IF (SM2%SINGLE_SPEC_INDEX==H2O_INDEX) THEN
      NL_DATA = NLMBDMIE
   ELSEIF (SPECIES(SM2%SINGLE_SPEC_INDEX)%ISFUEL) THEN
      NL_DATA = NFUELK
   ELSE
      NL_DATA = 1
   ENDIF
ENDIF

!     ALLOCATE ARRAYS

ALLOCATE(QSCA(1:NRDMIE,1:NLMBDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','QSCA',IZERO)
ALLOCATE(QABS(1:NRDMIE,1:NLMBDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','QABS',IZERO)
ALLOCATE(CHI_F(1:NRDMIE,1:NLMBDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','CHI_F',IZERO)
!
ALLOCATE(LMBDMIE(1:NLMBDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','LMBDMIE',IZERO)
ALLOCATE(LMBDWGHT(1:NLMBDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','LMBDWGHT',IZERO)
ALLOCATE(REAL_REF_INDX(1:NLMBDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','REAL_REF_INDX',IZERO)
ALLOCATE(CMPLX_REF_INDX(1:NLMBDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','CMPLX_REF_INDX',IZERO)

! RADIATIVE PROPERTIES

LMBDMIE(1:NLMBDMIE) = CPLXREF_WATER(1:NLMBDMIE,1)

IF (PRESENT(PARTICLE_CLASS)) THEN
   IF (LPC%RADIATIVE_PROPERTY_INDEX > 0) THEN
      DO NX = 1,NLMBDMIE
         CALL INTERPOLATE1D(TA%TABLE_DATA(:,1)*1.0E-6_EB,TA%TABLE_DATA(:,2),LMBDMIE(NX),REAL_REF_INDX(NX))
         CALL INTERPOLATE1D(TA%TABLE_DATA(:,1)*1.0E-6_EB,TA%TABLE_DATA(:,3),LMBDMIE(NX),CMPLX_REF_INDX(NX))
      ENDDO
   ELSEIF (LPC%Y_INDEX==H2O_INDEX) THEN
      REAL_REF_INDX              = 1.33_EB
      CMPLX_REF_INDX(1:NLMBDMIE) = CPLXREF_WATER(1:NLMBDMIE,2)
   ELSEIF (LPC%FUEL) THEN
      DO NX = 1,NLMBDMIE
         CALL INTERPOLATE1D(CPLXREF_FUEL(:,1)*1.0E-6_EB,CPLXREF_FUEL(:,2),LMBDMIE(NX),REAL_REF_INDX(NX))
         CALL INTERPOLATE1D(CPLXREF_FUEL(:,1)*1.0E-6_EB,CPLXREF_FUEL(:,3),LMBDMIE(NX),CMPLX_REF_INDX(NX))
      ENDDO
   ELSE
      REAL_REF_INDX              = LPC%REAL_REFRACTIVE_INDEX
      CMPLX_REF_INDX             = LPC%COMPLEX_REFRACTIVE_INDEX
   ENDIF
ELSE
   IF (SM2%SINGLE_SPEC_INDEX==H2O_INDEX) THEN
      REAL_REF_INDX              = 1.33_EB
      CMPLX_REF_INDX(1:NLMBDMIE) = CPLXREF_WATER(1:NLMBDMIE,2)
   ELSEIF (SPECIES(SM2%SINGLE_SPEC_INDEX)%ISFUEL) THEN
      DO NX = 1,NLMBDMIE
         CALL INTERPOLATE1D(CPLXREF_FUEL(:,1)*1.0E-6_EB,CPLXREF_FUEL(:,2),LMBDMIE(NX),REAL_REF_INDX(NX))
         CALL INTERPOLATE1D(CPLXREF_FUEL(:,1)*1.0E-6_EB,CPLXREF_FUEL(:,3),LMBDMIE(NX),CMPLX_REF_INDX(NX))
      ENDDO
   ELSE
      REAL_REF_INDX              = SPECIES(SM2%SINGLE_SPEC_INDEX)%REAL_REFRACTIVE_INDEX
      CMPLX_REF_INDX             = SPECIES(SM2%SINGLE_SPEC_INDEX)%COMPLEX_REFRACTIVE_INDEX
   ENDIF
ENDIF

CALL MIE_SCATTERING

!     GENERATE INTEGRATION WEIGHTS FOR LAMBDA

IF (NLMBDMIE == 1) THEN
   LMBDWGHT(1) = 1._EB
ELSE
   LMBDWGHT(1) = 0.5_EB*(LMBDMIE(2)-LMBDMIE(1))
   DO I = 2, NLMBDMIE-1
      LMBDWGHT(I) = 0.5_EB*(LMBDMIE(I+1) - LMBDMIE(I-1))
   ENDDO
   LMBDWGHT(NLMBDMIE) = 0.5_EB*(LMBDMIE(NLMBDMIE)-LMBDMIE(NLMBDMIE-1))
ENDIF

ALLOCATE(WQABS(0:MIE_NDG,1:NUMBER_SPECTRAL_BANDS))
CALL ChkMemErr('IRAD','WQABS',IZERO)
WQABS = 0._EB
ALLOCATE(WQSCA(0:MIE_NDG,1:NUMBER_SPECTRAL_BANDS))
CALL ChkMemErr('IRAD','WQSCA',IZERO)
WQSCA = 0._EB
ALLOCATE(R50(0:MIE_NDG))
CALL ChkMemErr('IRAD','R50',IZERO)
R50 = 0._EB

!     LOOP OVER ALL RADIATION BANDS

BANDLOOP: DO IBND = 1,NSB
   !
   IF (NSB == 1) THEN
      NLAMBDALOW = 1
      NLAMBDAHIGH = NLMBDMIE
   ELSE
      NLAMBDALOW  = MINLOC(LMBDMIE,MASK=LMBDMIE>=WL_LOW(IBND)*1.E-6)
      NLAMBDAHIGH = MAXLOC(LMBDMIE,MASK=LMBDMIE<=WL_HIGH(IBND)*1.0E-6)
   ENDIF

   !     LOOP OVER ALL PARTICLE SIZE GROUPS


   DRGROUPLOOP: DO ND = 1, MIE_NDG

      R50(ND) = EXP(DGROUP_A*REAL(ND,EB) + DGROUP_B)

      !     LOOP OVER WAVELENGTHS

      IBSUM = 0._EB

      DO J = NLAMBDALOW(1),NLAMBDAHIGH(1)
         IB = PLANCK(RADTMP, LMBDMIE(J)*1.0E6_EB)
         IBSUM = IBSUM + LMBDWGHT(J)*IB

         ASUM = 0._EB
         BSUM = 0._EB

         ! PROPERTIES AT D32

         CALL INTERPOLATE1D(RDMIE,QSCA(:,J),R50(ND),AVAL)
         CALL INTERPOLATE1D(RDMIE,CHI_F(:,J),R50(ND),BVAL)
         BVAL = (1._EB-BVAL)
         ASUM = AVAL*BVAL
         CALL INTERPOLATE1D(RDMIE,QABS(:,J),R50(ND),BVAL)
         BSUM = BVAL

         WQSCA(ND,IBND) = WQSCA(ND,IBND) + ASUM*LMBDWGHT(J)*IB
         WQABS(ND,IBND) = WQABS(ND,IBND) + BSUM*LMBDWGHT(J)*IB
      ENDDO

      !     NORMALIZE WITH BLACKBODY RADIATION

      WQSCA(ND,IBND)  = WQSCA(ND,IBND)/IBSUM
      WQABS(ND,IBND)  = WQABS(ND,IBND)/IBSUM

ENDDO DRGROUPLOOP
ENDDO BANDLOOP

IF (PRESENT(PARTICLE_CLASS)) THEN
   LPC%R50 = R50
   LPC%WQSCA = WQSCA
   LPC%WQABS = WQABS
ELSE
   SM2%R50 = R50
   SM2%WQSCA = WQSCA
   SM2%WQABS = WQABS
ENDIF

DEALLOCATE(WQABS)
DEALLOCATE(WQSCA)
DEALLOCATE(R50)

DEALLOCATE(RDMIE)
DEALLOCATE(QSCA)
DEALLOCATE(QABS)
DEALLOCATE(CHI_F)
DEALLOCATE(LMBDMIE)
DEALLOCATE(LMBDWGHT)
DEALLOCATE(REAL_REF_INDX)
DEALLOCATE(CMPLX_REF_INDX)

END SUBROUTINE MEAN_CROSS_SECTIONS



SUBROUTINE MIE_SCATTERING
!
!     CALCULATES THE SCATTERING AND ABSORPTION CROSS SECTIONS
!     AND CALCULATES FORWARD SCATTERING FRACTION BY INTEGRATING
!     THE SCATTERING PHASE FUNCTION.
!
! ----------------------------------------------------------------------
! -----------  SPECIFICATIONS FOR SUBROUTINE  MIEV0  ---------------
! ----------------------------------------------------------------------
INTEGER   MOMDIM
PARAMETER  ( MOMDIM = 200)
LOGICAL   ANYANG, PERFCT, PRNT( 2 )
INTEGER   IPOLZN, NMOM
REAL(EB)  GQSC, MIMCUT, PMOM( 0:MOMDIM, 4 ), SPIKE, QE, QS
REAL(EB), ALLOCATABLE :: XXX(:)
COMPLEX(EB) SFORW, SBACK, TFORW( 2 ), TBACK( 2 ), CREFIN
COMPLEX(EB), ALLOCATABLE :: S1(:), S2(:)

! --------------- LOCAL VARIABLES --------------------------------------

INTEGER   I, J, K, NX, IZERO
INTEGER   NLAMBDA, NRA, NMIEANG2
REAL(EB)  STMP, AIJ, FTMP, XX_MAX
REAL(EB)  MUMIN1, MUMIN2,THETALIM1, THETALIM2, MUDLOC
REAL(EB)  MU1, MU2, NU1, NU2, MUD0LOC, MUDPILOC, MUD1, MUD2, DMUD
REAL(EB),ALLOCATABLE :: XMU1(:), XNU1(:),XMU2(:), &
          ANGLE1(:), ANGLE2(:), MUD(:), MUDX(:),PWGHT(:), &
          PHSFUN(:), PFOR(:,:), MUD0(:,:), MUDPI(:,:)
! ----------------------------------------------------------------------

NRA    = NUMBER_RADIATION_ANGLES

!     MIEV-CODE VARIABLES

MIMCUT = 1.E-6_EB
PERFCT = .FALSE.
ANYANG = .TRUE.
!      IPOLZN = +1234
IPOLZN = 0
NMOM   = 0
PRNT   = .FALSE.

!     LIMIT FOR XX

XX_MAX = 15000.0_EB

!     INTEGRATION LIMITS

THETALIM1 = ACOS(1._EB - 2._EB/REAL(NRA))
MUMIN1 = COS(THETALIM1)
MUMIN2 = MUMIN1**2-(1._EB-MUMIN1**2)
THETALIM2 = ACOS(MUMIN2)
NMIEANG2 = NMIEANG*2

!     ALLOCATE LOCAL ARRAYS

ALLOCATE(XXX(1:NRDMIE),STAT=IZERO)
CALL CHKMEMERR('INIT','XXX',IZERO)
ALLOCATE(S1(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','S1',IZERO)
ALLOCATE(S2(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','S2',IZERO)
ALLOCATE(XMU1(1:NMIEANG),STAT=IZERO)
CALL CHKMEMERR('INIT','XMU1',IZERO)
ALLOCATE(XNU1(1:NMIEANG),STAT=IZERO)
CALL CHKMEMERR('INIT','XNU1',IZERO)
ALLOCATE(XMU2(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','XMU2',IZERO)
ALLOCATE(MUD(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','MUD',IZERO)
ALLOCATE(MUDX(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','MUDX',IZERO)
ALLOCATE(PWGHT(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','PWGHT',IZERO)
ALLOCATE(ANGLE1(1:NMIEANG),STAT=IZERO)
CALL CHKMEMERR('INIT','ANGLE1',IZERO)
ALLOCATE(ANGLE2(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','ANGLE2',IZERO)
ALLOCATE(PHSFUN(1:NMIEANG2),STAT=IZERO)
CALL CHKMEMERR('INIT','PHSFUN',IZERO)
ALLOCATE(PFOR(1:NMIEANG,1:NMIEANG),STAT=IZERO)
CALL CHKMEMERR('INIT','PFOR',IZERO)
ALLOCATE(MUD0(1:NMIEANG,1:NMIEANG),STAT=IZERO)
CALL CHKMEMERR('INIT','MUD0',IZERO)
ALLOCATE(MUDPI(1:NMIEANG,1:NMIEANG),STAT=IZERO)
CALL CHKMEMERR('INIT','MUDPI',IZERO)

!     CREATE SOLID ANGLE INTEGRATION ARRAYS

DO I = 1,NMIEANG
   ANGLE1(I) = THETALIM1*REAL(I-1)/REAL(NMIEANG-1)
   ANGLE1(I) = ANGLE1(I)*(1._EB-0.99*REAL(NMIEANG-I)/REAL(NMIEANG))
ENDDO
DO I = 1,NMIEANG
   XMU1(NMIEANG-I+1) = COS(ANGLE1(I))
   XNU1(NMIEANG-I+1) = SIN(ANGLE1(I))
ENDDO

!    CREATE PHASE FUNCTION INGTEGRATION ARRAYS

DO I = 1,NMIEANG2
   ANGLE2(I) = THETALIM2*REAL(I-1)/REAL(NMIEANG2-1)
   ANGLE2(I) = ANGLE2(I)*(1._EB-0.99*REAL(NMIEANG2-I)/REAL(NMIEANG2))
ENDDO
DO I = 1,NMIEANG2
   XMU2(NMIEANG2-I+1) = COS(ANGLE2(I))
ENDDO

!    CALCULATE PHASE FUNCTION INGETRATION LIMITS

DO J = 1,NMIEANG
   DO I = 1,NMIEANG
      MUD0(I,J)  = XMU1(I)*XMU1(J) + XNU1(I)*XNU1(J)
      MUDPI(I,J) = XMU1(I)*XMU1(J) - XNU1(I)*XNU1(J)
   ENDDO
ENDDO

!     CALCULATE PHASE FUNCTION INTEGRATION WEIGHTS

MU1 = 0.7_EB
MU2 = 0.9_EB
NU1 = SQRT(1-MU1**2)
NU2 = SQRT(1-MU2**2)
MUD0LOC  = MU1*MU2 + NU1*NU2
MUDPILOC = MU1*MU2 - NU1*NU2
MUD1 = MUDPILOC
MUD2 = MUD0LOC
DMUD = (MUD2-MUD1)/(NMIEANG2-1)
DO I = 1,NMIEANG2
   MUD(I) = MUD1+REAL(I-1)*DMUD
ENDDO
MUD(1)       = MUD(1)       + 0.25_EB*DMUD !EMPIRICAL
MUD(NMIEANG2) = MUD(NMIEANG2) - 0.25_EB*DMUD !EMPIRICAL
DO I = 1,NMIEANG2
   MUDX(I) = (MUD(I)-MUD1)/(MUD2-MUD1)
ENDDO
DO I = 1,NMIEANG2
   PWGHT(I) = DMUD/SQRT((NU1*NU2)**2-(MUD(I)-MU1*MU2)**2)
ENDDO
PWGHT(2)         = 0.5_EB*PWGHT(2)
PWGHT(NMIEANG2-1) = 0.5_EB*PWGHT(NMIEANG2-1)

!     LOOP OVER WAVELENGTH

LAMBDALOOP: DO NLAMBDA = 1, NLMBDMIE
!
   CREFIN = CMPLX( REAL_REF_INDX(NLAMBDA), CMPLX_REF_INDX(NLAMBDA), EB )

! CHOOSE PERFECTLY REFLECTING SPHERE, IF LARGE REAL INDEX IS GIVEN.

   IF (REAL_REF_INDX(NLAMBDA) > 10._EB) PERFCT = .TRUE.

!     LOOP OVER PARTICLE RADIUS

   RADIUSLOOP: DO NX = 1, NRDMIE

      XXX(NX) = MIN(XX_MAX,2._EB*PI*RDMIE(NX)/LMBDMIE(NLAMBDA))
      CALL MIEV0( XXX(NX), CREFIN, PERFCT, MIMCUT, ANYANG,   &
                  NMIEANG2, XMU2, NMOM, IPOLZN, MOMDIM, PRNT,  &
                  QE, QS, GQSC, &
                  PMOM, SFORW, SBACK, S1, S2, TFORW, TBACK, SPIKE )

      QSCA(NX,NLAMBDA) = QS
      QABS(NX,NLAMBDA) = QE-QS

!     CALCULATE SINGLE DROP PHASE FUNCTION

      IF (ABS(QS)>TWO_EPSILON_EB) THEN
         DO I = 1,NMIEANG2
            PHSFUN(I) = 2._EB*(ABS(S1(I))**2 + ABS(S2(I))**2 )
         ENDDO
         PHSFUN = PHSFUN/(QS*XXX(NX)**2)
      ELSE
         PHSFUN = 1.0_EB
      ENDIF

!     CALCULATE THE INNERMOST INTEGRAL OF THE FORWARD SCATTERING FRACTION

      PFOR = 0._EB
      DO J = 1,NMIEANG
         DO I = J,NMIEANG
            IF (ABS(MUD0(I,J)-MUDPI(I,J))<=TWO_EPSILON_EB) THEN
               CALL INTERPOLATE1D(XMU2,PHSFUN,MUD0(I,J),FTMP)
               PFOR(I,J) = PI*FTMP
            ELSE
               MUD1 = MUDPI(I,J)
               MUD2 = MUD0(I,J)
               STMP = 0.0_EB
               DO K = 1,NMIEANG2
                  MUDLOC = MUD1+MUDX(K)*(MUD2-MUD1)
                  CALL INTERPOLATE1D(XMU2,PHSFUN,MUDLOC,FTMP)
                  STMP = STMP + PWGHT(K)*FTMP
               ENDDO
               PFOR(I,J) = STMP
            ENDIF
         ENDDO
      ENDDO

!     CALCULATE THE TWO OUTER INTEGRALS OF THE FORWARD FRACTION

      STMP = 0._EB
      DO J = 1,NMIEANG-1
         DO I = J+1,NMIEANG-1
            AIJ = (XMU1(I+1)-XMU1(I))*(XMU1(J+1)-XMU1(J))/2._EB
            STMP = STMP +   (2._EB*PFOR(I,J)+PFOR(I+1,J)+PFOR(I,J+1)+2._EB*PFOR(I+1,J+1))*AIJ/3._EB
         ENDDO
      ENDDO
      DO I = 1,NMIEANG-1
         AIJ = ((XMU1(I+1)-XMU1(I))**2)/2._EB
         STMP = STMP + (PFOR(I,I)+PFOR(I+1,I)+PFOR(I+1,I+1))*AIJ/3._EB
      ENDDO
      CHI_F(NX,NLAMBDA) = 2._EB*STMP/(4._EB*PI/NRA)
   ENDDO RADIUSLOOP
ENDDO LAMBDALOOP

DEALLOCATE(XXX)
DEALLOCATE(S1)
DEALLOCATE(S2)
DEALLOCATE(XMU1)
DEALLOCATE(XNU1)
DEALLOCATE(XMU2)
DEALLOCATE(MUD)
DEALLOCATE(MUDX)
DEALLOCATE(PWGHT)
DEALLOCATE(ANGLE1)
DEALLOCATE(ANGLE2)
DEALLOCATE(PHSFUN)
DEALLOCATE(PFOR)
DEALLOCATE(MUD0)
DEALLOCATE(MUDPI)

END SUBROUTINE MIE_SCATTERING

SUBROUTINE MIEV0( XX, CREFIN, PERFCT, MIMCUT, ANYANG, NUMANG, XMU, &
                  NMOM, IPOLZN, MOMDIM, PRNT, QEXT, QSCA, GQSC, &
                  PMOM, SFORW, SBACK, S1, S2, TFORW, TBACK, &
                  SPIKE )

!     MIE SCATTERING FOR A SINGLE PARTICLE AND WAVELENGTH.
!     AUTHOR:  DR. WARREN J. WISCOMBE (WISCOMBE@CLIMATE.GSFC.NASA.GOV)
!         NASA GODDARD SPACE FLIGHT CENTER
!         CODE 913
!         GREENBELT, MD 20771

!     REFERENCES
!     ----------
!
!     (1) WISCOMBE, W., 1979: MIE SCATTERING CALCULATIONS--ADVANCES
!         IN TECHNIQUE AND FAST, VECTOR-SPEED COMPUTER CODES,
!         NCAR TECH NOTE TN-140+STR, NATIONAL CENTER FOR
!         ATMOSPHERIC RESEARCH, BOULDER, COLORADO (OUT OF PRINT
!         BUT AN UPDATED ELECTRONIC VERSION AVAILABLE)
!
!     (2) WISCOMBE, W., 1980: IMPROVED MIE SCATTERING ALGORITHMS,
!         APPL. OPT. 19, 1505-1509
!

!    COMPUTES MIE SCATTERING AND EXTINCTION EFFICIENCIES; ASYMMETRY
!    FACTOR;  FORWARD- AND BACKSCATTER AMPLITUDE;  SCATTERING
!    AMPLITUDES VS. SCATTERING ANGLE FOR INCIDENT POLARIZATION PARALLEL
!    AND PERPENDICULAR TO THE PLANE OF SCATTERING;
!    COEFFICIENTS IN THE LEGENDRE POLYNOMIAL EXPANSIONS OF EITHER THE
!    UNPOLARIZED PHASE FUNCTION OR THE POLARIZED PHASE MATRIX;
!    SOME QUANTITIES NEEDED IN POLARIZED RADIATIVE TRANSFER;  AND
!    INFORMATION ABOUT WHETHER OR NOT A RESONANCE HAS BEEN HIT.
!
!    INPUT AND OUTPUT VARIABLES ARE DESCRIBED IN FILE MIEV.DOC.
!    MANY STATEMENTS ARE ACCOMPANIED BY COMMENTS REFERRING TO
!    REFERENCES IN MIEV.DOC, NOTABLY THE NCAR MIE REPORT WHICH IS NOW
!    AVAILABLE ELECTRONICALLY AND WHICH IS REFERRED TO USING THE
!    SHORTHAND (RN), MEANING EQ. (N) OF THE REPORT.

!    CALLING TREE:
!
!        MIEV0
!            TESTMI
!                TSTBAD
!                MIPRNT
!                ERRMSG
!            CKINMI
!                WRTBAD
!                WRTDIM
!                ERRMSG
!            SMALL1
!            SMALL2
!            ERRMSG
!            BIGA
!                CONFRA
!                    ERRMSG
!            LPCOEF
!                LPCO1T
!                LPCO2T
!                ERRMSG
!            MIPRNT
!
!   I N P U T   V A R I A B L E S
!   -----------------------------
!
!  ( EVEN IF AN INPUT VARIABLE IS NOT NEEDED FOR A PARTICULAR
!    APPLICATION, MAKE SURE IT HAS A LEGITIMATE VALUE THAT CAN
!    BE WRITTEN OUT AND READ IN -- NO INDEFINITES, ETC. )
!
!  XX        MIE SIZE PARAMETER ( 2 * PI * RADIUS / WAVELENGTH )
!
!  CREFIN    COMPLEX REFRACTIVE INDEX ( IMAG PART CAN BE + OR -,
!            BUT INTERNALLY A NEGATIVE IMAGINARY INDEX IS ASSUMED ).
!            IF IMAG PART IS - ,  SCATTERING AMPLITUDES AS IN VAN
!            DE HULST ARE RETURNED;  IF IMAG PART IS + , COMPLEX
!            CONJUGATES OF THOSE SCATTERING AMPLITUDES ARE RETURNED
!            (THE LATTER IS THE CONVENTION IN PHYSICS).
!            ** NOTE ** IN THE 'PERFECT' CASE, SCATTERING AMPLITUDES
!            IN THE VAN DE HULST (REF. 6 ABOVE) CONVENTION WILL
!            AUTOMATICALLY BE RETURNED UNLESS  IM(CREFIN)  IS
!            POSITIVE;  OTHERWISE, CREFIN PLAYS NO ROLE.
!
!  PERFCT    TRUE, ASSUME REFRACTIVE INDEX IS INFINITE AND USE
!            SPECIAL CASE FORMULAS FOR MIE COEFFICIENTS  'A'
!            AND  'B'  ( SEE KERKER, M., THE SCATTERING OF
!            LIGHT AND OTHER ELECTROMAGNETIC RADIATION, P. 90 ).
!            THIS IS SOMETIMES CALLED THE 'TOTALLY REFLECTING',
!            SOMETIMES THE 'PERFECTLY CONDUCTING' CASE.
!            ( SEE CREFIN FOR ADDITIONAL INFORMATION )
!
!  MIMCUT    (POSITIVE) VALUE BELOW WHICH IMAGINARY REFRACTIVE
!            INDEX IS REGARDED AS ZERO (COMPUTATION PROCEEDS
!            FASTER FOR ZERO IMAGINARY INDEX)
!
!  ANYANG    TRUE, ANY ANGLES WHATSOEVER MAY BE INPUT THROUGH
!            XMU.  FALSE, THE ANGLES ARE MONOTONE INCREASING
!            AND MIRROR SYMMETRIC ABOUT 90 DEGREES (THIS OPTION
!            IS ADVANTAGEOUS BECAUSE THE SCATTERING AMPLITUDES
!            S1,S2 FOR THE ANGLES BETWEEN 90 AND 180 DEGREES
!            ARE EVALUABLE FROM SYMMETRY RELATIONS, AND HENCE
!            ARE OBTAINED WITH LITTLE ADDED COMPUTATIONAL COST.)
!
!  NUMANG    NO. OF ANGLES AT WHICH SCATTERING AMPLITUDES
!            S1,S2 ARE TO BE EVALUATED  ( SET = 0 TO SKIP
!            CALCULATION OF S1,S2 ).  MAKE SURE NUMANG DOES
!            NOT EXCEED THE PARAMETER MAXANG IN THE PROGRAM.
!
!  XMU(N)    COSINES OF ANGLES ( N = 1 TO NUMANG ) AT WHICH S1,S2
!            ARE TO BE EVALUATED.  IF ANYANG = FALSE, THEN
!
!             (A) THE ANGLES MUST BE MONOTONE INCREASING AND
!                 MIRROR SYMMETRIC ABOUT 90 DEGREES (IF 90-A IS
!                 AN ANGLE, THEN 90+A MUST BE ALSO)
!
!             (B) IF NUMANG IS ODD, 90 DEGREES MUST BE AMONG
!                 THE ANGLES
!
!  NMOM       HIGHEST LEGENDRE MOMENT PMOM TO CALCULATE,
!             NUMBERING FROM ZERO ( NMOM = 0 PREVENTS
!             CALCULATION OF PMOM )
!
!  IPOLZN     POSITIVE, COMPUTE LEGENDRE MOMENTS PMOM FOR THE
!                       MUELLER MATRIX ELEMENTS DETERMINED BY THE
!                       DIGITS OF IPOLZN, WITH 1 REFERRING TO M1,
!                       2 TO M2, 3 TO S21, AND 4 TO D21 (REF. 3).
!                       E.G., IF IPOLZN = 14 THEN ONLY MOMENTS FOR
!                       M1 AND D21 WILL BE RETURNED.
!
!             0,        COMPUTE LEGENDRE MOMENTS PMOM FOR THE
!                       UNPOLARIZED UNNORMALIZED PHASE FUNCTION.
!
!             NEGATIVE, COMPUTE LEGENDRE MOMENTS PMOM FOR THE
!                       SEKERA PHASE QUANTITIES DETERMINED BY THE
!                       DIGITS OF ABS(IPOLZN), WITH 1 REFERRING TO
!                       R1, 2 TO R2, 3 TO R3, AND 4 TO R4 (REF. 4).
!                       E.G., IF IPOLZN = -14 THEN ONLY MOMENTS FOR
!                       R1 AND R4 WILL BE RETURNED.
!
!             ( NOT USED IF  NMOM = 0 )
!
!  MOMDIM     DETERMINES FIRST DIMENSION OF PMOM, WHICH IS DIMENSIONED
!             INTERNALLY AS PMOM( 0:MOMDIM, * ) (SECOND DIMENSION MUST
!             BE THE LARGER OF UNITY AND THE HIGHEST DIGIT IN
!             IPOLZN; IF NOT, SERIOUS ERRORS WILL OCCUR).
!             MUST BE GIVEN A VALUE, EVEN IF  NMOM = 0.  MINIMUM: 1.
!
!  PRT(L)     PRINT FLAGS (LOGICAL).  L = 1  PRINTS  S1,S2, THEIR
!             SQUARED ABSOLUTE VALUES, AND DEGREE OF POLARIZATION,
!             PROVIDED NUMANG IS NON-ZERO.   L = 2  PRINTS ALL
!             OUTPUT VARIABLES OTHER THAN  S1,S2.
!
!
! O U T P U T   V A R I A B L E S
! -------------------------------
!
!  QEXT      (REAL) EXTINCTION EFFICIENCY FACTOR  ( REF. 2, EQ. 1A )
!
!  QSCA      (REAL) SCATTERING EFFICIENCY FACTOR  ( REF. 2, EQ. 1B )
!
!  GQSC      (REAL) ASYMMETRY FACTOR TIMES SCATTERING EFFICIENCY
!            ( REF. 2, EQ. 1C )  ( ALLOWS CALCULATION OF RADIATION
!            PRESSURE EFFICIENCY FACTOR  QPR = QEXT - GQSC )
!
!  =====================================================================
!  ==== NOTE --  S1, S2, SFORW, SBACK, TFORW, AND TBACK ARE CALCULATED
!  ====          INTERNALLY FOR NEGATIVE IMAGINARY REFRACTIVE INDEX;
!  ====          FOR POSITIVE IMAGINARY INDEX, THEIR COMPLEX CONJUGATES
!  ====          ARE TAKEN BEFORE THEY ARE RETURNED, TO CORRESPOND TO
!  ====          CUSTOMARY USAGE IN SOME PARTS OF PHYSICS ( IN PARTI-
!  ====          CULAR, IN PAPERS ON CAM APPROXIMATIONS TO MIE THEORY ).
!  =====================================================================
!
!  S1(N),    (COMPLEX) MIE SCATTERING AMPLITUDES AT ANGLES SPECIFIED
!  S2(N)     BY XMU(N) ( N=1 TO NUMANG )  ( REF. 2, EQS. 1D-E ).
!
!  SFORW     (COMPLEX) FORWARD-SCATTERING AMPLITUDE S1 AT
!            0 DEGREES.  ( S2(0 DEG) = S1(0 DEG) )
!
!  SBACK     (COMPLEX) BACKSCATTERING AMPLITUDE S1 AT
!            180 DEGREES.   ( S2(180 DEG) = - S1(180 DEG) )
!
!  TFORW(I)  (COMPLEX) VALUES OF
!
!                I=1:  T1 = ( S2 - (MU)*S1 ) / ( 1 - MU**2 )
!                I=2:  T2 = ( S1 - (MU)*S2 ) / ( 1 - MU**2 )
!
!            AT ANGLE THETA = 0 ( MU = COS(THETA) = 1 ), WHERE THE
!            EXPRESSIONS ON THE RIGHT-HAND SIDE ARE INDETERMINATE.
!            ( THESE QUANTITIES ARE REQUIRED FOR DOING POLARIZED
!            RADIATIVE TRANSFER (REF. 4, APPENDIX). )
!  TBACK(I)  (COMPLEX) VALUES OF  T1 (FOR I=1) OR  T2 (FOR I=2) AT
!            ANGLE  THETA = 180 DEGREES ( MU = COS(THETA) = - 1 ).
!
!  SPIKE     (REAL) MAGNITUDE OF THE SMALLEST DENOMINATOR OF
!            EITHER MIE COEFFICIENT (A-SUB-N OR B-SUB-N),
!            TAKEN OVER ALL TERMS IN THE MIE SERIES PAST
!            N = SIZE PARAMETER XX.  VALUES OF SPIKE BELOW
!            ABOUT 0.3 SIGNIFY A RIPPLE SPIKE, SINCE THESE
!            SPIKES ARE PRODUCED BY ABNORMALLY SMALL DENOMINATORS
!            IN THE MIE COEFFICIENTS (NORMAL DENOMINATORS ARE OF
!            ORDER UNITY OR HIGHER).  DEFAULTS TO 1.0 WHEN NOT
!            ON A SPIKE.  DOES NOT IDENTIFY ALL RESONANCES
!            (WE ARE STILL WORKING ON THAT).
!
! PMOM(M,NP) (REAL) MOMENTS  M = 0 TO NMOM  OF UNNORMALIZED NP-TH
!            PHASE QUANTITY PQ  ( MOMENTS WITH  M > 2*NTRM  ARE
!            ZERO, WHERE  NTRM = NO. TERMS IN MIE SERIES =
!            XX + 4*XX**1/3 + 1 ) :
!
!              PQ( MU, NP ) = SUM( M=0 TO INFINITY ) ( (2M+1)
!                                * PMOM( M,NP ) * P-SUB-M( MU ) )
!
!            WHERE  MU = COS( SCATTERING ANGLE )
!                   P-SUB-M = M-TH LEGENDRE POLYNOMIAL
!
!            AND THE DEFINITION OF 'PQ' IS AS FOLLOWS:
!
!            IPOLZN>0:  PQ(MU,1) = CABS( S1(MU) )**2
!                          PQ(MU,2) = CABS( S2(MU) )**2
!                          PQ(MU,3) = RE( S1(MU)*CONJG( S2(MU) ) )
!                          PQ(MU,4) = - IM( S1(MU)*CONJG( S2(MU) ) )
!                          ( CALLED M1, M2, S21, D21 IN LITERATURE )
!
!            IPOLZN=0:  PQ(MU,1) = ( CABS(S1)**2 + CABS(S2)**2 ) / 2
!                       ( THE UNNORMALIZED PHASE FUNCTION )
!
!            IPOLZN<0:  PQ(MU,1) = CABS( T1(MU) )**2
!                          PQ(MU,2) = CABS( T2(MU) )**2
!                          PQ(MU,3) = RE( T1(MU)*CONJG( T2(MU) ) )
!                          PQ(MU,4) = - IM( T1(MU)*CONJG( T2(MU) ) )
!                          ( CALLED R1, R2, R3, R4 IN LITERATURE )
!
!            THE SIGN OF THE 4TH PHASE QUANTITY IS A SOURCE OF
!            CONFUSION.  IT FLIPS IF THE COMPLEX CONJUGATES OF
!            S1,S2  OR  T1,T2  ARE USED, AS OCCURS WHEN A
!            REFRACTIVE INDEX WITH POSITIVE IMAGINARY PART IS
!            USED (SEE DISCUSSION BELOW).  THE DEFINITION ABOVE
!            IS CONSISTENT WITH A NEGATIVE IMAGINARY PART.
!
!            ** WARNING **  MAKE SURE THE SECOND DIMENSION OF PMOM
!            IN THE CALLING PROGRAM IS AT LEAST AS LARGE AS THE
!            ABSOLUTE VALUE OF IPOLZN.
!
!            FOR SMALL ENOUGH VALUES OF XX, OR LARGE ENOUGH VALUES
!            OF M,  PMOM  WILL TEND TO UNDERFLOW.  THUS, IT IS
!            UNWISE TO ASSUME THE VALUES RETURNED ARE NON-ZERO AND,
!            FOR EXAMPLE, TO DIVIDE SOME QUANTITY BY THEM.
!


!      I N T E R N A L   V A R I A B L E S
!      -----------------------------------

!  AN,BN           MIE COEFFICIENTS A-SUB-N, B-SUB-N ( REF. 1, EQ. 16 )
!  ANM1,BNM1       MIE COEFFICIENTS  A-SUB-(N-1),
!                     B-SUB-(N-1);  USED IN GQSC SUM
!  ANP             COEFFS. IN S+ EXPANSION ( REF. 2, P. 1507 )
!  BNP             COEFFS. IN S- EXPANSION ( REF. 2, P. 1507 )
!  ANPM            COEFFS. IN S+ EXPANSION ( REF. 2, P. 1507 )
!                     WHEN  MU  IS REPLACED BY  - MU
!  BNPM            COEFFS. IN S- EXPANSION ( REF. 2, P. 1507 )
!                     WHEN  MU  IS REPLACED BY  - MU
!  CALCMO(K)       TRUE, CALCULATE MOMENTS FOR K-TH PHASE QUANTITY
!                     (DERIVED FROM IPOLZN)
!  CBIGA(N)        BESSEL FUNCTION RATIO A-SUB-N (REF. 2, EQ. 2)
!                     ( COMPLEX VERSION )
!  CDENAN,         (COMPLEX) DENOMINATORS OF AN,BN
!   CDENBN
!  CIOR            COMPLEX INDEX OF REFRACTION WITH NEGATIVE
!                     IMAGINARY PART (VAN DE HULST CONVENTION)
!  CIORIV          1 / CIOR
!  COEFF           ( 2N + 1 ) / ( N ( N + 1 ) )
!  CSUM1,2         TEMPORARY SUM VARIABLES FOR TFORW, TBACK
!  FN              FLOATING POINT VERSION OF LOOP INDEX FOR
!                     MIE SERIES SUMMATION
!  LITA,LITB(N)    MIE COEFFICIENTS AN, BN, SAVED IN ARRAYS FOR
!                     USE IN CALCULATING LEGENDRE MOMENTS PMOM
!  MAXTRM          MAX. POSSIBLE NO. OF TERMS IN MIE SERIES
!  MM              (-1)^(N+1), WHERE N IS MIE SERIES SUM INDEX
!  MIM             MAGNITUDE OF IMAGINARY REFRACTIVE INDEX
!  MRE             REAL PART OF REFRACTIVE INDEX
!  MAXANG          MAX. POSSIBLE VALUE OF INPUT VARIABLE NUMANG
!  NANGD2          (NUMANG+1)/2 ( NO. OF ANGLES IN 0-90 DEG; ANYANG=F )
!  NOABS           TRUE, SPHERE NON-ABSORBING (DETERMINED BY MIMCUT)
!  NP1DN           ( N + 1 ) / N
!  NPQUAN          HIGHEST-NUMBERED PHASE QUANTITY FOR WHICH MOMENTS ARE
!                     TO BE CALCULATED (THE LARGEST DIGIT IN IPOLZN
!                     IF  IPOLZN /= 0)
!  NTRM            NO. OF TERMS IN MIE SERIES
!  PASS1           TRUE ON FIRST ENTRY, FALSE THEREAFTER; FOR SELF-TEST
!  PIN(J)          ANGULAR FUNCTION PI-SUB-N ( REF. 2, EQ. 3 )
!                     AT J-TH ANGLE
!  PINM1(J)        PI-SUB-(N-1) ( SEE PIN ) AT J-TH ANGLE
!  PSINM1          RICATTI-BESSEL FUNCTION PSI-SUB-(N-1), ARGUMENT XX
!  PSIN            RICATTI-BESSEL FUNCTION PSI-SUB-N OF ARGUMENT XX
!                     ( REF. 1, P. 11 FF. )
!  RBIGA(N)        BESSEL FUNCTION RATIO A-SUB-N (REF. 2, EQ. 2)
!                     ( REAL VERSION, FOR WHEN IMAG REFRAC INDEX = 0 )
!  RIORIV          1 / MRE
!  RN              1 / N
!  RTMP            (REAL) TEMPORARY VARIABLE
!  SP(J)           S+  FOR J-TH ANGLE  ( REF. 2, P. 1507 )
!  SM(J)           S-  FOR J-TH ANGLE  ( REF. 2, P. 1507 )
!  SPS(J)          S+  FOR (NUMANG+1-J)-TH ANGLE ( ANYANG=FALSE )
!  SMS(J)          S-  FOR (NUMANG+1-J)-TH ANGLE ( ANYANG=FALSE )
!  TAUN            ANGULAR FUNCTION TAU-SUB-N ( REF. 2, EQ. 4 )
!                     AT J-TH ANGLE
!  TCOEF           N ( N+1 ) ( 2N+1 ) (FOR SUMMING TFORW,TBACK SERIES)
!  TWONP1          2N + 1
!  YESANG          TRUE IF SCATTERING AMPLITUDES ARE TO BE CALCULATED
!  ZETNM1          RICATTI-BESSEL FUNCTION  ZETA-SUB-(N-1) OF ARGUMENT
!                     XX  ( REF. 2, EQ. 17 )
!  ZETN            RICATTI-BESSEL FUNCTION  ZETA-SUB-N OF ARGUMENT XX
! ----------------------------------------------------------------------
!
!
! ----------------------------------------------------------------------
! --------  I / O SPECIFICATIONS FOR SUBROUTINE MIEV0  -----------------
! ----------------------------------------------------------------------
LOGICAL     ANYANG, PERFCT, PRNT(*)
INTEGER     IPOLZN, MOMDIM, NUMANG, NMOM
REAL(EB)    GQSC, MIMCUT, PMOM( 0:MOMDIM, * ), QEXT, QSCA, SPIKE,XMU(*), XX
COMPLEX(EB) CREFIN, SFORW, SBACK, S1(*), S2(*), TFORW(*), TBACK(*)
! ----------------------------------------------------------------------
!
!                                  ** NOTE --  MAXTRM = 10100  IS NECES-
!                                  ** SARY TO DO SOME OF THE TEST PROBS,
!                                  ** BUT 1100 IS SUFFICIENT FOR MOST
!                                  ** CONCEIVABLE APPLICATIONS
!     .. PARAMETERS ..

INTEGER     MAXANG
PARAMETER   ( MAXANG = 5000 )
INTEGER     MAXTRM
!      PARAMETER ( MAXTRM = 10100 )
PARAMETER   ( MAXTRM = 16000 ) ! WORKS FOR FDS
REAL(EB)    ONETHR
PARAMETER   ( ONETHR = 1._EB / 3._EB )

!     .. LOCAL SCALARS ..

LOGICAL     NOABS, PASS1, YESANG, GT100
INTEGER     I, J, N, NANGD2, NPQUAN, NTRM
REAL(EB)    CHIN, CHINM1, COEFF, DENAN, DENBN, FN, MIM, MM, MRE, &
            NP1DN, PSIN, PSINM1, RATIO, RIORIV, RN, RTMP, TAUN, &
            TCOEF, TWONP1, XINV
COMPLEX(EB) AN, ANM1, ANP, ANPM, BN, BNM1, BNP, BNPM, CDENAN,&
            CDENBN, CIOR, CIORIV, CSUM1, CSUM2, ZET, ZETN, ZETNM1

!     .. LOCAL ARRAYS ..

LOGICAL     CALCMO( 4 )
REAL(EB), ALLOCATABLE :: PIN(:), PINM1(:)
REAL(EB) RBIGA( MAXTRM )
COMPLEX(EB) CBIGA( MAXTRM ), LITA( MAXTRM ), LITB( MAXTRM )
COMPLEX(EB), ALLOCATABLE :: SM(:), SMS(:), SP(:), SPS(:)

SAVE      PASS1
DATA      PASS1 / .TRUE. /
!
ALLOCATE(PIN(1:NUMANG))
ALLOCATE(PINM1(1:NUMANG))
ALLOCATE(SM(1:NUMANG))
ALLOCATE(SMS(1:(NUMANG+1)/2))
ALLOCATE(SP(1:NUMANG))
ALLOCATE(SPS(1:(NUMANG+1)/2))

!                    ** SAVE SOME INPUT VARIABLES AND REPLACE THEM
!                    ** WITH VALUES NEEDED TO DO THE SELF-TEST

IF( PASS1 ) CALL TESTMI( .FALSE., XX, CREFIN, MIMCUT, PERFCT, &
                         ANYANG, NMOM, IPOLZN, NUMANG, XMU, QEXT, &
                         QSCA, GQSC, SFORW, SBACK, S1, S2, TFORW, &
                         TBACK, PMOM, MOMDIM )
MAIN_MIEV: DO
!                                        ** CHECK INPUT AND CALCULATE
!                                        ** CERTAIN VARIABLES FROM INPUT
CALL CKINMI( NUMANG, MAXANG, XX, PERFCT, CREFIN, MOMDIM, NMOM, &
             IPOLZN, ANYANG, XMU, CALCMO, NPQUAN )

IF( PERFCT .AND. XX<=0.1_EB ) THEN
!                                            ** USE TOTALLY-REFLECTING
!                                            ** SMALL-PARTICLE LIMIT
   CALL SMALL1( XX, NUMANG, XMU, QEXT, QSCA, GQSC, SFORW, SBACK, &
                  S1, S2, TFORW, TBACK, LITA, LITB )
   NTRM = 2
ELSE
   NOABS = .TRUE.
   GT100 = .FALSE.

   IF( .NOT.PERFCT ) THEN
      CIOR = CREFIN
      IF( AIMAG(CIOR)>0.0_EB ) CIOR = CONJG( CIOR )
      MRE    = REAL( CIOR )
      MIM    = -AIMAG( CIOR )
      NOABS  = MIM<=MIMCUT
      CIORIV = 1.0_EB / CIOR
      RIORIV = 1.0_EB / MRE

      IF( XX*MAX( 1._EB, ABS(CIOR) )<=0.1_EB ) THEN
!                                    ** USE GENERAL-REFRACTIVE-INDEX
!                                    ** SMALL-PARTICLE LIMIT
         CALL SMALL2( XX, CIOR, MIM>MIMCUT, NUMANG, XMU, QEXT, &
                           QSCA, GQSC, SFORW, SBACK, S1, S2, TFORW, &
                           TBACK, LITA, LITB )
         NTRM = 2
         GT100 = .TRUE.
      END IF
   END IF
   GT100IF: IF (.NOT. GT100) THEN

      NANGD2 = ( NUMANG + 1 ) / 2
      YESANG = NUMANG>0
   !                             ** NUMBER OF TERMS IN MIE SERIES; EQ R50
      IF( XX<=8.0_EB ) THEN
      NTRM = INT(XX + 4._EB*XX**ONETHR + 1._EB)
      ELSE IF( XX<4200._EB ) THEN
      NTRM = INT(XX + 4.05*XX**ONETHR + 2._EB)
      ELSE
      NTRM =INT( XX + 4._EB*XX**ONETHR + 2._EB)
      END IF
      IF( NTRM+1 > MAXTRM )CALL ERRMSG('MIEV0--PARAMETER MAXTRM TOO SMALL',.TRUE.)
   !                            ** CALCULATE LOGARITHMIC DERIVATIVES OF
   !                            ** J-BESSEL-FCN., A-SUB-(1 TO NTRM)
      IF( .NOT.PERFCT ) CALL BIGA( CIOR, XX, NTRM, NOABS, YESANG, RBIGA, CBIGA )
   !                            ** INITIALIZE RICATTI-BESSEL FUNCTIONS
   !                            ** (PSI,CHI,ZETA)-SUB-(0,1) FOR UPWARD
   !                            ** RECURRENCE ( EQ. R19 )
      XINV   = 1.0_EB / XX
      PSINM1 = SIN( XX )
      CHINM1 = COS( XX )
      PSIN   = PSINM1*XINV - CHINM1
      CHIN   = CHINM1*XINV + PSINM1
      ZETNM1 = CMPLX( PSINM1, CHINM1, EB )
      ZETN   = CMPLX( PSIN, CHIN, EB )
   !                                     ** INITIALIZE PREVIOUS COEFFI-
   !                                     ** CIENTS FOR GQSC SERIES
      ANM1 = ( 0.0_EB, 0.0_EB )
      BNM1 = ( 0.0_EB, 0.0_EB )
   !                             ** INITIALIZE ANGULAR FUNCTION  PI
   !                             ** AND SUMS FOR S+, S- ( REF. 2, P. 1507 )
      IF( ANYANG ) THEN
      DO J = 1, NUMANG
   !                             ** EQ. R39
            PINM1( J ) = 0.0_EB
            PIN( J ) = 1.0_EB

               SP( J ) = ( 0.0_EB, 0.0_EB )
               SM( J ) = ( 0.0_EB, 0.0_EB )
         END DO
      ELSE
         DO J = 1, NANGD2
   !                          ** EQ. R39
            PINM1( J ) = 0.0_EB
            PIN( J ) = 1.0_EB
            SP( J ) = ( 0.0_EB, 0.0_EB )
            SM( J ) = ( 0.0_EB, 0.0_EB )
            SPS( J ) = ( 0.0_EB, 0.0_EB )
            SMS( J ) = ( 0.0_EB, 0.0_EB )
         END DO
      END IF
   !                       ** INITIALIZE MIE SUMS FOR EFFICIENCIES, ETC.
      QSCA  = 0.0_EB
      GQSC  = 0.0_EB
      SFORW = ( 0._EB, 0._EB )
      SBACK = ( 0._EB, 0._EB )
      CSUM1 = ( 0._EB, 0._EB )
      CSUM2 = ( 0._EB, 0._EB )
   !
   ! ---------  LOOP TO SUM MIE SERIES  -----------------------------------
      MM     = +1.0_EB
      SPIKE  = 1.0_EB
      DO N = 1, NTRM
   !                           ** COMPUTE VARIOUS NUMERICAL COEFFICIENTS
         FN     = N
         RN     = 1.0_EB / FN
         NP1DN  = 1.0_EB + RN
         TWONP1 = 2*N + 1
         COEFF  = TWONP1 / ( FN * ( N + 1 ) )
         TCOEF  = TWONP1 * ( FN * ( N + 1 ) )
   !                          ** CALCULATE MIE SERIES COEFFICIENTS
         IF( PERFCT ) THEN
   !                                 ** TOTALLY-REFLECTING CASE; EQ R/A.1,2
            AN = ( ( FN*XINV )*PSIN - PSINM1 ) / ( ( FN*XINV )*ZETN - ZETNM1 )
            BN = PSIN / ZETN
         ELSE IF( NOABS ) THEN
   !                                      ** NO-ABSORPTION CASE; EQ (R16)
            CDENAN = ( RIORIV*RBIGA(N) + ( FN*XINV ) ) * ZETN - ZETNM1
            AN   = ( ( RIORIV*RBIGA(N) + ( FN*XINV ) ) * PSIN - PSINM1 ) / CDENAN
            CDENBN = ( MRE*RBIGA(N) + ( FN*XINV ) ) * ZETN - ZETNM1
            BN   = ( ( MRE*RBIGA(N) + ( FN*XINV ) ) * PSIN - PSINM1 ) / CDENBN
         ELSE
   !                                       ** ABSORPTIVE CASE; EQ (R16)
            CDENAN = ( CIORIV*CBIGA( N ) + ( FN*XINV ) )*ZETN - ZETNM1
            CDENBN =   ( CIOR*CBIGA( N ) + ( FN*XINV ) )*ZETN - ZETNM1
            AN   = ( ( CIORIV*CBIGA( N ) + ( FN*XINV ) )*PSIN - PSINM1 ) / CDENAN
            BN     = ( ( CIOR*CBIGA( N ) + ( FN*XINV ) )*PSIN - PSINM1 ) / CDENBN
   !                                         ** EQ (R7)
            QSCA   = QSCA + TWONP1*( SQ( AN ) + SQ( BN ) )
         END IF
   !                       ** SAVE MIE COEFFICIENTS FOR PMOM CALCULATION
         LITA( N ) = AN
         LITB( N ) = BN
   !
         IF( .NOT.PERFCT .AND. N>XX ) THEN
   !                                               ** FLAG RESONANCE SPIKES
            DENAN  = ABS( CDENAN )
            DENBN  = ABS( CDENBN )
   !                                                   ** EQ. R/B.9
            RATIO  = DENAN / DENBN
   !                                                   ** EQ. R/B.10
            IF( RATIO<=0.2_EB .OR. RATIO>=5.0_EB ) SPIKE = MIN( SPIKE, DENAN, DENBN )
         END IF
   !                                  ** INCREMENT MIE SUMS FOR NON-ANGLE-
   !                                  ** DEPENDENT QUANTITIES
   !                                                   ** EQ. R/B.2
         SFORW = SFORW + TWONP1*( AN + BN )
   !                                                   ** EQ. R/B.5,6
         CSUM1 = CSUM1 + TCOEF *( AN - BN )
   !                                                   ** EQ. R/B.1
         SBACK = SBACK + ( MM*TWONP1 )*( AN - BN )
   !                                                   ** EQ. R/B.7,8
         CSUM2 = CSUM2 + ( MM*TCOEF ) *( AN + BN )
   !                                         ** EQ (R8)
         GQSC  = GQSC  + (FN - RN) * REAL(ANM1 * CONJG(AN) + BNM1 * CONJG(BN) ) + COEFF * REAL(AN * CONJG(BN))
         IF( YESANG ) THEN
   !                                      ** PUT MIE COEFFICIENTS IN FORM
   !                                      ** NEEDED FOR COMPUTING S+, S-
   !                                      ** ( EQ R10 )
         ANP = COEFF*( AN + BN )
         BNP = COEFF*( AN - BN )
   !                                      ** INCREMENT MIE SUMS FOR S+, S-
   !                                      ** WHILE UPWARD RECURSING
   !                                      ** ANGULAR FUNCTIONS PI AND TAU
         IF( ANYANG ) THEN
   !                                         ** ARBITRARY ANGLES
   !                                              ** VECTORIZABLE LOOP
            DO J = 1, NUMANG
   !                                                 ** EQ. (R37B)
               RTMP = ( XMU(J) * PIN(J) ) - PINM1( J )
   !                                                 ** EQ. (R38B)
               TAUN   = FN * RTMP - PINM1( J )
   !                                                   ** EQ (R10)
               SP( J ) = SP( J ) + ANP * ( PIN( J ) + TAUN )
               SM( J ) = SM( J ) + BNP * ( PIN( J ) - TAUN )
               PINM1( J ) = PIN( J )
   !                                                 ** EQ. R37C
               PIN( J ) = ( XMU( J ) * PIN( J ) ) + NP1DN * RTMP
            END DO
         ELSE
   !                                  ** ANGLES SYMMETRIC ABOUT 90 DEGREES
            ANPM = MM*ANP
            BNPM = MM*BNP
   !                                          ** VECTORIZABLE LOOP
            DO J = 1, NANGD2
   !                                                 ** EQ. (R37B)
               RTMP = ( XMU(J) * PIN(J) ) - PINM1( J )
   !                                                 ** EQ. (R38B)
               TAUN = FN * RTMP - PINM1( J )
   !                                                 ** EQ (R10,12)
               SP ( J ) = SP ( J ) + ANP * ( PIN( J ) + TAUN )
               SMS( J ) = SMS( J ) + BNPM *( PIN( J ) + TAUN )
               SM ( J ) = SM ( J ) + BNP * ( PIN( J ) - TAUN )
               SPS( J ) = SPS( J ) + ANPM *( PIN( J ) - TAUN )
               PINM1( J ) = PIN( J )
         !                                                 ** EQ. R37C
               PIN( J ) = ( XMU(J) * PIN(J) ) + NP1DN * RTMP
            END DO
         END IF
         END IF
   !                          ** UPDATE RELEVANT QUANTITIES FOR NEXT
   !                          ** PASS THROUGH LOOP
         MM   = - MM
         ANM1 = AN
         BNM1 = BN
   !                           ** UPWARD RECURRENCE FOR RICATTI-BESSEL
   !                           ** FUNCTIONS ( EQ. R17 )
         ZET    = ( TWONP1*XINV ) * ZETN - ZETNM1
         ZETNM1 = ZETN
         ZETN   = ZET
         PSINM1 = PSIN
         PSIN   = REAL( ZETN )
      END DO
   ! ---------- END LOOP TO SUM MIE SERIES --------------------------------
   !
   !                                         ** EQ (R6)
      QEXT = 2._EB / XX**2*REAL( SFORW )
      IF( PERFCT .OR. NOABS ) THEN
         QSCA = QEXT
      ELSE
         QSCA = 2._EB/ XX**2 * QSCA
      END IF
      GQSC   = 4._EB/ XX**2 * GQSC
      SFORW  = 0.5_EB*SFORW
      SBACK  = 0.5_EB*SBACK
      TFORW( 1 ) =  0.5_EB*SFORW - 0.125_EB*CSUM1
      TFORW( 2 ) =  0.5_EB*SFORW + 0.125_EB*CSUM1
      TBACK( 1 ) = -0.5_EB*SBACK + 0.125_EB*CSUM2
      TBACK( 2 ) =  0.5_EB*SBACK + 0.125_EB*CSUM2
      IF( YESANG ) THEN
   !                                ** RECOVER SCATTERING AMPLITUDES
   !                                ** FROM S+, S- ( EQ (R11) )
         IF( ANYANG ) THEN
   !                                         ** VECTORIZABLE LOOP
            DO J = 1, NUMANG
   !                                                  ** EQ (R11)
               S1( J ) = 0.5_EB*( SP( J ) + SM( J ) )
               S2( J ) = 0.5_EB*( SP( J ) - SM( J ) )
            END DO
         ELSE
   !                                        ** VECTORIZABLE LOOP
            DO J = 1, NANGD2
   !                                                   ** EQ (R11)
               S1( J ) = 0.5_EB*( SP( J ) + SM( J ) )
               S2( J ) = 0.5_EB*( SP( J ) - SM( J ) )
            END DO
   !                                        ** VECTORIZABLE LOOP
            DO  J = 1, NANGD2
               S1( NUMANG + 1 - J ) = 0.5_EB*( SPS( J ) + SMS( J ) )
               S2( NUMANG + 1 - J ) = 0.5_EB*( SPS( J ) - SMS( J ) )
            END DO
         END IF
      END IF
   END IF GT100IF
   !                                 ** CALCULATE LEGENDRE MOMENTS
END IF
IF( NMOM>0 ) CALL LPCOEF( NTRM, NMOM, IPOLZN, MOMDIM, CALCMO,NPQUAN, LITA, LITB, PMOM )
IF( AIMAG( CREFIN )>0.0_EB ) THEN
!                                         ** TAKE COMPLEX CONJUGATES
!                                         ** OF SCATTERING AMPLITUDES
   SFORW = CONJG( SFORW )
   SBACK = CONJG( SBACK )
   DO I = 1, 2
      TFORW( I ) = CONJG( TFORW( I ) )
      TBACK( I ) = CONJG( TBACK( I ) )
   END DO
   DO J = 1, NUMANG
      S1( J ) = CONJG( S1( J ) )
      S2( J ) = CONJG( S2( J ) )
   END DO
END IF

IF( PASS1 ) THEN
!                           ** COMPARE TEST CASE RESULTS WITH
!                           ** CORRECT ANSWERS AND ABORT IF BAD;
!                           ** OTHERWISE RESTORE USER INPUT AND PROCEED
   CALL TESTMI( .TRUE., XX, CREFIN, MIMCUT, PERFCT, ANYANG, NMOM, &
                  IPOLZN, NUMANG, XMU, QEXT, QSCA, GQSC, SFORW, &
                  SBACK, S1, S2, TFORW, TBACK, PMOM, MOMDIM )
   PASS1  = .FALSE.
   CYCLE MAIN_MIEV
END IF
EXIT MAIN_MIEV
END DO MAIN_MIEV

IF( PRNT( 1 ) .OR. PRNT( 2 ) )  &
  CALL MIPRNT( PRNT, XX, PERFCT, CREFIN, NUMANG, XMU, QEXT, &
               QSCA, GQSC, NMOM, IPOLZN, MOMDIM, CALCMO, PMOM, &
               SFORW, SBACK, TFORW, TBACK, S1, S2 )
!
DEALLOCATE(PIN)
DEALLOCATE(PINM1)
DEALLOCATE(SM)
DEALLOCATE(SMS)
DEALLOCATE(SP)
DEALLOCATE(SPS)

RETURN

END SUBROUTINE MIEV0

SUBROUTINE CKINMI( NUMANG, MAXANG, XX, PERFCT, CREFIN, MOMDIM, NMOM, IPOLZN, ANYANG, XMU, CALCMO, NPQUAN )

!        CHECK FOR BAD INPUT TO MIEV0 AND CALCULATE CALCMO, NPQUAN

!     ROUTINES CALLED :  ERRMSG, WRTBAD, WRTDIM

!     .. SCALAR ARGUMENTS ..

LOGICAL     ANYANG, PERFCT
INTEGER     IPOLZN, MAXANG, MOMDIM, NMOM, NPQUAN, NUMANG
REAL(EB)    XX
COMPLEX(EB) CREFIN
!     ..
!     .. ARRAY ARGUMENTS ..

LOGICAL     CALCMO( * )
REAL(EB)    XMU( * )
!     ..
!     .. LOCAL SCALARS ..

CHARACTER(4) :: STRING
LOGICAL     INPERR
INTEGER     I, IP, J, L

INPERR = .FALSE.

IF( NUMANG>MAXANG ) INPERR = WRTDIM( 'MAXANG', NUMANG )
IF( NUMANG<0 ) INPERR = WRTBAD( 'NUMANG' )

IF( XX<0._EB ) INPERR = WRTBAD( 'XX' )

IF( .NOT.PERFCT .AND. REAL( CREFIN )<=0._EB )INPERR = WRTBAD( 'CREFIN' )

IF( MOMDIM<0 ) INPERR = WRTBAD( 'MOMDIM' )


IF( NMOM/=0 ) THEN

   IF( NMOM<0 .OR. NMOM>MOMDIM ) INPERR = WRTBAD( 'NMOM' )
   IF( ABS( IPOLZN )>4444 ) INPERR = WRTBAD( 'IPOLZN' )

   NPQUAN = 0

   DO L = 1, 4
      CALCMO( L ) = .FALSE.
   END DO

   IF( IPOLZN/=0 ) THEN
   !                                 ** PARSE OUT IPOLZN INTO ITS DIGITS
   !                                 ** TO FIND WHICH PHASE QUANTITIES ARE
   !                                 ** TO HAVE THEIR MOMENTS CALCULATED
      WRITE( STRING, '(I4)' ) ABS( IPOLZN )
      DO J = 1, 4
         IP = ICHAR( STRING( J:J ) ) - ICHAR( '0' )
         IF( IP>=1 .AND. IP<=4 ) CALCMO( IP ) = .TRUE.
         IF( IP==0 .OR. ( IP>=5 .AND. IP<=9 ) ) INPERR = WRTBAD( 'IPOLZN' )
         NPQUAN = MAX( NPQUAN, IP )
      END DO
   END IF
END IF

IF( ANYANG ) THEN
!                                ** ALLOW FOR SLIGHT IMPERFECTIONS IN
!                                ** COMPUTATION OF COSINE
   DO I = 1, NUMANG
      IF( XMU( I )<-1.00001_EB .OR. XMU( I )>1.00001_EB ) INPERR = WRTBAD( 'XMU' )
   END DO
   ELSE
   DO I = 1, ( NUMANG + 1 ) / 2
      IF( XMU( I )<-0.00001_EB .OR. XMU( I )>1.00001_EB ) INPERR = WRTBAD( 'XMU' )
   END DO
END IF

IF( INPERR ) CALL ERRMSG( 'MIEV0--INPUT ERROR(S).  ABORTING...', .TRUE. )
IF( XX>20000.0_EB .OR. REAL( CREFIN )>10.0_EB .OR. ABS( AIMAG( CREFIN ) )>10.0_EB ) &
    CALL ERRMSG( 'MIEV0--XX OR CREFIN OUTSIDE TESTED RANGE',.FALSE.)
RETURN
END SUBROUTINE CKINMI

SUBROUTINE LPCOEF( NTRM, NMOM, IPOLZN, MOMDIM, CALCMO, NPQUAN, A, B, PMOM )
!         CALCULATE LEGENDRE POLYNOMIAL EXPANSION COEFFICIENTS (ALSO
!         CALLED MOMENTS) FOR PHASE QUANTITIES ( REF. 5 FORMULATION )

!     INPUT:  NTRM                    NUMBER TERMS IN MIE SERIES
!             NMOM, IPOLZN, MOMDIM    MIEV0 ARGUMENTS
!             CALCMO                  FLAGS CALCULATED FROM IPOLZN
!             NPQUAN                  DEFINED IN MIEV0
!             A, B                    MIE SERIES COEFFICIENTS
!
!     OUTPUT: PMOM                   LEGENDRE MOMENTS (MIEV0 ARGUMENT)
!
!     ROUTINES CALLED :  ERRMSG, LPCO1T, LPCO2T
!
!     *** NOTES ***
!
!         (1)  EQS. 2-5 ARE IN ERROR IN DAVE, APPL. OPT. 9,
!         1888 (1970).  EQ. 2 REFERS TO M1, NOT M2;  EQ. 3 REFERS TO
!         M2, NOT M1.  IN EQS. 4 AND 5, THE SUBSCRIPTS ON THE SECOND
!         TERM IN SQUARE BRACKETS SHOULD BE INTERCHANGED.
!
!         (2)  THE GENERAL-CASE LOGIC IN THIS SUBROUTINE WORKS CORRECTLY
!         IN THE TWO-TERM MIE SERIES CASE, BUT SUBROUTINE LPCO2T
!         IS CALLED INSTEAD, FOR SPEED.
!
!         (3)  SUBROUTINE  LPCO1T, TO DO THE ONE-TERM CASE, IS NEVER
!         CALLED WITHIN THE CONTEXT OF MIEV0, BUT IS INCLUDED FOR
!         COMPLETE GENERALITY.
!
!         (4)  SOME IMPROVEMENT IN SPEED IS OBTAINABLE BY COMBINING THE
!         310- AND 410-LOOPS, IF MOMENTS FOR BOTH THE THIRD AND FOURTH
!         PHASE QUANTITIES ARE DESIRED, BECAUSE THE THIRD PHASE QUANTITY
!         IS THE REAL PART OF A COMPLEX SERIES, WHILE THE FOURTH PHASE
!         QUANTITY IS THE IMAGINARY PART OF THAT VERY SAME SERIES.  BUT
!         MOST USERS ARE NOT INTERESTED IN THE FOURTH PHASE QUANTITY,
!         WHICH IS RELATED TO CIRCULAR POLARIZATION, SO THE PRESENT
!         SCHEME IS USUALLY MORE EFFICIENT.
!
!
!           ** DEFINITIONS OF LOCAL VARIABLES ***
!      AM(M)       NUMERICAL COEFFICIENTS  A-SUB-M-SUPER-L
!                     IN DAVE, EQS. 1-15, AS SIMPLIFIED IN REF. 5.
!      BI(I)       NUMERICAL COEFFICIENTS  B-SUB-I-SUPER-L
!                     IN DAVE, EQS. 1-15, AS SIMPLIFIED IN REF. 5.
!      BIDEL(I)    1/2 BI(I) TIMES FACTOR CAPITAL-DEL IN DAVE
!      C,D()       ARRAYS C AND D IN DAVE, EQS. 16-17 (MUELLER FORM),
!                     CALCULATED USING RECURRENCE DERIVED IN REF. 5
!      CS,DS()     ARRAYS C AND D IN REF. 4, EQS. A5-A6 (SEKERA FORM),
!                     CALCULATED USING RECURRENCE DERIVED IN REF. 5
!      EVENL       TRUE FOR EVEN-NUMBERED MOMENTS;  FALSE OTHERWISE
!      IDEL        1 + LITTLE-DEL  IN DAVE
!      MAXTRM      MAX. NO. OF TERMS IN MIE SERIES
!      MAXMOM      MAX. NO. OF NON-ZERO MOMENTS
!      NUMMOM      NUMBER OF NON-ZERO MOMENTS
!      RECIP(K)    1 / K

!     .. PARAMETERS ..

INTEGER     MAXTRM, MAXMOM, MXMOM2, MAXRCP
PARAMETER ( MAXTRM=1001, MAXMOM = 2*MAXTRM, MXMOM2 = MAXMOM / 2, MAXRCP = 4*MAXTRM + 2 )
!     ..
!     .. SCALAR ARGUMENTS ..

INTEGER   IPOLZN, MOMDIM, NMOM, NPQUAN, NTRM
!     ..
!     .. ARRAY ARGUMENTS ..

LOGICAL     CALCMO( * )
REAL(EB)    PMOM( 0:MOMDIM, * )
COMPLEX(EB) A( * ), B( * )
!     ..
!     .. LOCAL SCALARS ..

LOGICAL     EVENL, PASS1
INTEGER     I, IDEL, IMAX, J, K, L, LD2, M, MMAX, NUMMOM
REAL(EB)    SUM
!     ..
!     .. LOCAL ARRAYS ..

REAL(EB)    AM( 0:MAXTRM ), BI( 0:MXMOM2 ), BIDEL( 0:MXMOM2 ), RECIP( MAXRCP )
COMPLEX(EB) C( MAXTRM ), CS( MAXTRM ), D( MAXTRM ), DS( MAXTRM )
SAVE      PASS1, RECIP
DATA      PASS1 / .TRUE. /

IF( PASS1 ) THEN
   DO  K = 1, MAXRCP
      RECIP( K ) = 1.0_EB / K
   END DO
   PASS1  = .FALSE.
END IF

DO J = 1, MAX( 1, NPQUAN )
   DO L = 0, NMOM
      PMOM( L, J ) = 0.0_EB
   END DO
END DO

IF( NTRM==1 ) THEN
   CALL LPCO1T( NMOM, IPOLZN, MOMDIM, CALCMO, A, B, PMOM )
   RETURN
ELSE IF( NTRM==2 ) THEN
   CALL LPCO2T( NMOM, IPOLZN, MOMDIM, CALCMO, A, B, PMOM )
   RETURN
END IF

IF( NTRM + 2>MAXTRM )CALL ERRMSG('LPCOEF--PARAMETER MAXTRM TOO SMALL',.TRUE.)
!                                     ** CALCULATE MUELLER C, D ARRAYS
C( NTRM + 2 ) = ( 0._EB, 0._EB )
D( NTRM + 2 ) = ( 0._EB, 0._EB )
C( NTRM + 1 ) = ( 1._EB - RECIP( NTRM+1 ) ) * B( NTRM )
D( NTRM + 1 ) = ( 1._EB - RECIP( NTRM+1 ) ) * A( NTRM )
C( NTRM ) = ( RECIP( NTRM ) + RECIP( NTRM+1 ) ) * A( NTRM ) + ( 1._EB - RECIP( NTRM ) )*B( NTRM-1 )
D( NTRM ) = ( RECIP( NTRM ) + RECIP( NTRM+1 ) ) * B( NTRM ) + ( 1._EB - RECIP( NTRM ) )*A( NTRM-1 )
DO K = NTRM-1, 2, -1
   C( K ) = C( K+2 ) - ( 1._EB + RECIP(K+1) ) * B( K+1 )+ ( RECIP(K) + RECIP(K+1) ) * A( K ) &
                        + ( 1._EB - RECIP(K) ) * B( K-1 )
   D( K ) = D( K+2 ) - ( 1._EB + RECIP(K+1) ) * A( K+1 )  + ( RECIP(K) + RECIP(K+1) ) * B( K ) &
                        + ( 1._EB - RECIP(K) ) * A( K-1 )
END DO

C( 1 ) = C( 3 ) + 1.5_EB * ( A( 1 ) - B( 2 ) )
D( 1 ) = D( 3 ) + 1.5_EB * ( B( 1 ) - A( 2 ) )

IF( IPOLZN>=0 ) THEN
   DO K = 1, NTRM + 2
      C( K ) = ( 2*K - 1 ) * C( K )
      D( K ) = ( 2*K - 1 ) * D( K )
   END DO
ELSE
!                                    ** COMPUTE SEKERA C AND D ARRAYS
   CS( NTRM + 2 ) = ( 0._EB, 0._EB )
   DS( NTRM + 2 ) = ( 0._EB, 0._EB )
   CS( NTRM + 1 ) = ( 0._EB, 0._EB )
   DS( NTRM + 1 ) = ( 0._EB, 0._EB )

   DO K = NTRM, 1, -1
      CS( K ) = CS( K+2 ) + ( 2*K + 1 ) * ( C( K+1 ) - B( K ) )
      DS( K ) = DS( K+2 ) + ( 2*K + 1 ) * ( D( K+1 ) - A( K ) )
   END DO

   DO K = 1, NTRM + 2
      C( K ) = ( 2*K - 1 ) * CS( K )
      D( K ) = ( 2*K - 1 ) * DS( K )
   END DO

END IF

IF( IPOLZN<0 ) NUMMOM = MIN( NMOM, 2*NTRM - 2 )
IF( IPOLZN>=0 ) NUMMOM = MIN( NMOM, 2*NTRM )

IF( NUMMOM>MAXMOM )  CALL ERRMSG('LPCOEF--PARAMETER MAXTRM TOO SMALL',.TRUE.)
!
!                          ** LOOP OVER MOMENTS
L240: DO L = 0, NUMMOM
   LD2 = L / 2
   EVENL  = MOD( L, 2 )==0
!                                    ** CALCULATE NUMERICAL COEFFICIENTS
!                                    ** A-SUB-M AND B-SUB-I IN DAVE
!                                    ** DOUBLE-SUMS FOR MOMENTS
   IF( L==0 ) THEN
      IDEL = 1
      DO M = 0, NTRM
         AM( M ) = 2.0_EB * RECIP( 2*M + 1 )
      END DO
      BI( 0 ) = 1.0_EB
   ELSE IF( EVENL ) THEN
      IDEL = 1
      DO M = LD2, NTRM
         AM( M ) = ( 1._EB + RECIP( 2*M - L + 1 ) ) * AM( M )
      END DO
      DO I = 0, LD2 - 1
         BI( I ) = ( 1._EB - RECIP( L - 2*I ) ) * BI( I )
      END DO
      BI( LD2 ) = ( 2._EB - RECIP( L ) ) * BI( LD2 - 1 )
   ELSE
   IDEL = 2
      DO M = LD2, NTRM
         AM( M ) = ( 1._EB - RECIP( 2*M + L + 2 ) ) * AM( M )
      END DO

      DO I = 0, LD2
         BI( I ) = ( 1._EB - RECIP( L + 2*I + 1 ) ) * BI( I )
      END DO
   END IF
!                                     ** ESTABLISH UPPER LIMITS FOR SUMS
!                                     ** AND INCORPORATE FACTOR CAPITAL-
!                                     ** DEL INTO B-SUB-I
   MMAX = NTRM - IDEL
   IF( IPOLZN>=0 ) MMAX = MMAX + 1
   IMAX = MIN( LD2, MMAX - LD2 )
   IF( IMAX<0 ) EXIT L240
   DO I = 0, IMAX
      BIDEL( I ) = BI( I )
   END DO
   IF( EVENL ) BIDEL( 0 ) = 0.5_EB*BIDEL( 0 )
!                                    ** PERFORM DOUBLE SUMS JUST FOR
!                                    ** PHASE QUANTITIES DESIRED BY USER
   IF( IPOLZN==0 ) THEN
      DO I = 0, IMAX
!                                           ** VECTORIZABLE LOOP
         SUM = 0.0_EB
         DO M = LD2, MMAX - I
            SUM = SUM + AM(M) * ( REAL(C(M-I+1) * CONJG(C(M+I+IDEL)))+ REAL(D(M-I+1) * CONJG(D(M+I+IDEL))))
         END DO
         PMOM( L, 1 ) = PMOM( L, 1 ) + BIDEL( I ) * SUM
      END DO
      PMOM( L, 1 ) = 0.5_EB*PMOM( L, 1 )
      CYCLE L240
   END IF

   IF( CALCMO( 1 ) ) THEN
      DO I = 0, IMAX
         SUM = 0.0_EB
!                                          ** VECTORIZABLE LOOP
         DO M = LD2, MMAX - I
            SUM = SUM + AM( M ) * REAL( C(M-I+1) * CONJG( C(M+I+IDEL) ) )
         END DO
         PMOM( L, 1 ) = PMOM( L, 1 ) + BIDEL( I ) * SUM
      END DO
   END IF

   IF( CALCMO( 2 ) ) THEN
      DO I = 0, IMAX
         SUM = 0.0_EB
!                                           ** VECTORIZABLE LOOP
         DO M = LD2, MMAX - I
            SUM = SUM + AM( M ) * REAL( D(M-I+1) * CONJG( D(M+I+IDEL) ) )
         END DO
         PMOM( L, 2 ) = PMOM( L, 2 ) + BIDEL( I ) * SUM
      END DO
   END IF

   IF( CALCMO( 3 ) ) THEN
      DO I = 0, IMAX
         SUM = 0.0_EB
!                                           ** VECTORIZABLE LOOP
         DO M = LD2, MMAX - I
            SUM = SUM + AM(M) *( REAL(C(M-I+1) * CONJG(D(M+I+IDEL))) + REAL(C(M+I+IDEL) * CONJG(D(M-I+1))))
         END DO
         PMOM( L, 3 ) = PMOM( L, 3 ) + BIDEL( I ) * SUM
      END DO
      PMOM( L, 3 ) = 0.5_EB*PMOM( L, 3 )
   END IF

   IF( CALCMO( 4 ) ) THEN
      DO I = 0, IMAX
         SUM= 0.0_EB
!                                         ** VECTORIZABLE LOOP
         DO M = LD2, MMAX - I
            SUM = SUM + AM(M) *(AIMAG(C(M-I+1) * CONJG(D(M+I+IDEL)))+ AIMAG(C(M+I+IDEL) * CONJG(D(M-I+1))))
         END DO
         PMOM( L, 4 ) = PMOM( L, 4 ) + BIDEL( I ) * SUM
      END DO

      PMOM( L, 4 ) = - 0.5_EB * PMOM( L, 4 )
   END IF

END DO L240

RETURN
END SUBROUTINE LPCOEF

SUBROUTINE LPCO1T( NMOM, IPOLZN, MOMDIM, CALCMO, A, B, PMOM )
!
!         CALCULATE LEGENDRE POLYNOMIAL EXPANSION COEFFICIENTS (ALSO
!         CALLED MOMENTS) FOR PHASE QUANTITIES IN SPECIAL CASE WHERE
!         NO. TERMS IN MIE SERIES = 1
!
!        INPUT:  NMOM, IPOLZN, MOMDIM     MIEV0 ARGUMENTS
!                CALCMO                   FLAGS CALCULATED FROM IPOLZN
!                A(1), B(1)               MIE SERIES COEFFICIENTS
!
!        OUTPUT: PMOM                     LEGENDRE MOMENTS

!     .. SCALAR ARGUMENTS ..

INTEGER     IPOLZN, MOMDIM, NMOM
!     ..
!     .. ARRAY ARGUMENTS ..

LOGICAL     CALCMO( * )
REAL(EB)    PMOM( 0:MOMDIM, * )
COMPLEX(EB) A( * ), B( * )
!     ..
!     .. LOCAL SCALARS ..

INTEGER     L, NUMMOM
REAL(EB)    A1SQ, B1SQ
COMPLEX(EB) A1B1C

A1SQ   = SQ( A( 1 ) )
B1SQ   = SQ( B( 1 ) )
A1B1C  = A( 1 ) * CONJG( B( 1 ) )


IF( IPOLZN<0 ) THEN
   IF( CALCMO( 1 ) ) PMOM( 0, 1 ) = 2.25_EB*B1SQ
   IF( CALCMO( 2 ) ) PMOM( 0, 2 ) = 2.25_EB*A1SQ
   IF( CALCMO( 3 ) ) PMOM( 0, 3 ) = 2.25_EB*REAL( A1B1C )
   IF( CALCMO( 4 ) ) PMOM( 0, 4 ) = 2.25_EB*AIMAG( A1B1C )
ELSE
   NUMMOM = MIN( NMOM, 2 )
   !                             ** LOOP OVER MOMENTS
   L10: DO L = 0, NUMMOM
      IF( IPOLZN==0 ) THEN
         IF( L==0 ) PMOM( L, 1 ) = 1.5_EB*( A1SQ + B1SQ )
         IF( L==1 ) PMOM( L, 1 ) = 1.5_EB*REAL( A1B1C )
         IF( L==2 ) PMOM( L, 1 ) = 0.15_EB*( A1SQ + B1SQ )
         CYCLE L10
      END IF

      IF( CALCMO( 1 ) ) THEN
         IF( L==0 ) PMOM( L, 1 ) = 2.25_EB*( A1SQ + B1SQ / 3.)
         IF( L==1 ) PMOM( L, 1 ) = 1.5_EB*REAL( A1B1C )
         IF( L==2 ) PMOM( L, 1 ) = 0.3_EB*B1SQ
      END IF

      IF( CALCMO( 2 ) ) THEN
         IF( L==0 ) PMOM( L, 2 ) = 2.25_EB*( B1SQ + A1SQ / 3. )
         IF( L==1 ) PMOM( L, 2 ) = 1.5_EB*REAL( A1B1C )
         IF( L==2 ) PMOM( L, 2 ) = 0.3_EB*A1SQ
      END IF

      IF( CALCMO( 3 ) ) THEN
         IF( L==0 ) PMOM( L, 3 ) = 3.0_EB*REAL( A1B1C )
         IF( L==1 ) PMOM( L, 3 ) = 0.75_EB*( A1SQ + B1SQ )
         IF( L==2 ) PMOM( L, 3 ) = 0.3_EB*REAL( A1B1C )
      END IF

      IF( CALCMO( 4 ) ) THEN
         IF( L==0 ) PMOM( L, 4 ) = -1.5_EB*AIMAG( A1B1C )
         IF( L==1 ) PMOM( L, 4 ) = 0.0_EB
         IF( L==2 ) PMOM( L, 4 ) = 0.3_EB*AIMAG( A1B1C )
      END IF

   END DO L10
END IF

RETURN
END SUBROUTINE LPCO1T

SUBROUTINE LPCO2T( NMOM, IPOLZN, MOMDIM, CALCMO, A, B, PMOM )
!         CALCULATE LEGENDRE POLYNOMIAL EXPANSION COEFFICIENTS (ALSO
!         CALLED MOMENTS) FOR PHASE QUANTITIES IN SPECIAL CASE WHERE
!         NO. TERMS IN MIE SERIES = 2

!        INPUT:  NMOM, IPOLZN, MOMDIM     MIEV0 ARGUMENTS
!                CALCMO                   FLAGS CALCULATED FROM IPOLZN
!                A(1-2), B(1-2)           MIE SERIES COEFFICIENTS
!
!        OUTPUT: PMOM                     LEGENDRE MOMENTS


!     .. SCALAR ARGUMENTS ..

INTEGER     IPOLZN, MOMDIM, NMOM
!     ..
!     .. ARRAY ARGUMENTS ..

LOGICAL     CALCMO( * )
REAL(EB)    PMOM( 0:MOMDIM, * )
COMPLEX(EB) A( * ), B( * )
!     ..
!     .. LOCAL SCALARS ..

INTEGER     L, NUMMOM
REAL(EB)    A2SQ, B2SQ, PM1, PM2
COMPLEX(EB) A2C, B2C, CA, CAC, CAT, CB, CBC, CBT, CG, CH

CA   = 3.*A( 1 ) - 5.*B( 2 )
CAT  = 3.*B( 1 ) - 5.*A( 2 )
CAC  = CONJG( CA )
A2SQ = SQ( A( 2 ) )
B2SQ = SQ( B( 2 ) )
A2C  = CONJG( A( 2 ) )
B2C  = CONJG( B( 2 ) )

IF( IPOLZN<0 ) THEN
!                                   ** LOOP OVER SEKERA MOMENTS
   NUMMOM = MIN( NMOM, 2 )
   DO L = 0, NUMMOM
      IF( CALCMO( 1 ) ) THEN
         IF( L==0 ) PMOM( L, 1 ) = 0.25_EB * ( SQ( CAT ) + (100._EB/3._EB)* B2SQ )
         IF( L==1 ) PMOM( L, 1 ) = (5._EB/3._EB)*REAL( CAT*B2C )
         IF( L==2 ) PMOM( L, 1 ) = (10._EB/3._EB)*B2SQ
      END IF

      IF( CALCMO( 2 ) ) THEN
         IF( L==0 ) PMOM( L, 2 ) = 0.25_EB * ( SQ( CA ) + (100._EB/3._EB) * A2SQ )
         IF( L==1 ) PMOM( L, 2 ) = (5._EB/3._EB)*REAL( CA*A2C )
         IF( L==2 ) PMOM( L, 2 ) = (10._EB/3._EB)*A2SQ
      END IF

      IF( CALCMO( 3 ) ) THEN
         IF( L==0 ) PMOM( L, 3 ) = 0.25_EB * REAL( CAT * CAC  + (100._EB/3._EB) * B(2) * A2C )
         IF( L==1 ) PMOM( L, 3 ) = 5._EB/6._EB* REAL( B(2)*CAC + CAT*A2C )
         IF( L==2 ) PMOM( L, 3 ) = 10._EB/3._EB* REAL( B(2)*A2C )
      END IF

      IF( CALCMO( 4 ) ) THEN
         IF( L==0 ) PMOM( L, 4 ) = -0.25_EB * AIMAG( CAT * CAC + (100._EB/3._EB)* B(2) * A2C )
         IF( L==1 ) PMOM( L, 4 ) = -5._EB/ 6._EB*  AIMAG( B(2)*CAC + CAT*A2C )
         IF( L==2 ) PMOM( L, 4 ) = -10._EB/ 3._EB* AIMAG( B(2)*A2C )
      END IF

   END DO

ELSE

   CB  = 3._EB*B( 1 ) + 5._EB*A( 2 )
   CBT = 3._EB*A( 1 ) + 5._EB*B( 2 )
   CBC = CONJG( CB )
   CG  = ( CBC*CBT + 10._EB*( CAC*A( 2 ) + B2C*CAT ) ) / 3._EB
   CH  = 2._EB*( CBC*A( 2 ) + B2C*CBT )

   !                               ** LOOP OVER MUELLER MOMENTS
   NUMMOM = MIN( NMOM, 4 )

   L20: DO L = 0, NUMMOM

      IF( IPOLZN==0 .OR. CALCMO( 1 ) ) THEN
         IF( L==0 ) PM1 = 0.25_EB*SQ( CA ) + SQ( CB ) / 12._EB + (5._EB/3._EB)*REAL( CA*B2C ) + 5._EB*B2SQ
         IF( L==1 ) PM1 = REAL( CB * ( CAC / 6._EB+ B2C ) )
         IF( L==2 ) PM1 = SQ( CB ) / 30._EB+ (20._EB/7._EB)*B2SQ + (2._EB/3._EB)*REAL( CA*B2C )
         IF( L==3 ) PM1 = (2._EB/7.) * REAL( CB*B2C )
         IF( L==4 ) PM1 = (40._EB/63._EB) * B2SQ
         IF( CALCMO( 1 ) ) PMOM( L, 1 ) = PM1
      END IF

      IF( IPOLZN==0 .OR. CALCMO( 2 ) ) THEN
         IF( L==0 ) PM2 = 0.25_EB*SQ( CAT ) + SQ( CBT ) / 12._EB + ( 5._EB/ 3._EB) * REAL( CAT*A2C ) + 5._EB*A2SQ
         IF( L==1 ) PM2 = REAL( CBT * ( CONJG( CAT ) / 6._EB+ A2C ) )
         IF( L==2 ) PM2 = SQ( CBT ) / 30._EB + ( 20._EB/7._EB) * A2SQ + ( 2._EB/3._EB) * REAL( CAT*A2C )
         IF( L==3 ) PM2 = (2._EB/7._EB) * REAL( CBT*A2C )
         IF( L==4 ) PM2 = (40._EB/63._EB) * A2SQ
         IF( CALCMO( 2 ) ) PMOM( L, 2 ) = PM2
      END IF

      IF( IPOLZN==0 ) THEN
         PMOM( L, 1 ) = 0.5_EB*( PM1 + PM2 )
         CYCLE L20
      END IF

      IF( CALCMO( 3 ) ) THEN
         IF( L==0 ) PMOM( L, 3 ) = 0.25_EB * REAL( CAC*CAT + CG + 20._EB* B2C * A(2) )
         IF( L==1 ) PMOM( L, 3 ) = REAL( CAC*CBT + CBC*CAT + 3._EB*CH ) / 12.
         IF( L==2 ) PMOM( L, 3 ) = 0.1_EB * REAL( CG + (200._EB/7._EB) * B2C * A(2) )
         IF( L==3 ) PMOM( L, 3 ) = REAL( CH ) / 14._EB
         IF( L==4 ) PMOM( L, 3 ) = 40._EB/63._EB* REAL( B2C*A(2) )
      END IF

      IF( CALCMO( 4 ) ) THEN
         IF( L==0 ) PMOM( L, 4 ) = 0.25_EB * AIMAG( CAC*CAT + CG + 20._EB* B2C * A(2) )
         IF( L==1 ) PMOM( L, 4 ) = AIMAG( CAC*CBT + CBC*CAT + 3._EB*CH ) / 12._EB
         IF( L==2 ) PMOM( L, 4 ) = 0.1_EB * AIMAG( CG + (200._EB/7._EB) * B2C * A(2) )
         IF( L==3 ) PMOM( L, 4 ) = AIMAG( CH ) / 14._EB
         IF( L==4 ) PMOM( L, 4 ) = 40._EB/63._EB* AIMAG( B2C*A(2) )
      END IF

   END DO L20

END IF

RETURN
END SUBROUTINE LPCO2T

SUBROUTINE BIGA( CIOR, XX, NTRM, NOABS, YESANG, RBIGA, CBIGA )

!        CALCULATE LOGARITHMIC DERIVATIVES OF J-BESSEL-FUNCTION
!
!     INPUT :  CIOR, XX, NTRM, NOABS, YESANG  (DEFINED IN MIEV0)
!
!    OUTPUT :  RBIGA OR CBIGA  (DEFINED IN MIEV0)
!
!    ROUTINES CALLED :  CONFRA
!
!
!    INTERNAL VARIABLES :
!
!       CONFRA     VALUE OF LENTZ CONTINUED FRACTION FOR CBIGA(NTRM),
!                     USED TO INITIALIZE DOWNWARD RECURRENCE
!
!       DOWN       = TRUE, USE DOWN-RECURRENCE.  FALSE, DO NOT.
!
!       F1,F2,F3   ARITHMETIC STATEMENT FUNCTIONS USED IN DETERMINING
!                     WHETHER TO USE UP-  OR DOWN-RECURRENCE
!                     ( REF. 2, EQS. 6-8 )
!
!       MRE        REAL REFRACTIVE INDEX
!       MIM        IMAGINARY REFRACTIVE INDEX
!
!       REZINV     1 / ( MRE * XX ); TEMPORARY VARIABLE FOR RECURRENCE
!       ZINV       1 / ( CIOR * XX ); TEMPORARY VARIABLE FOR RECURRENCE
!
!     .. SCALAR ARGUMENTS ..

LOGICAL     NOABS, YESANG
INTEGER     NTRM
REAL(EB)    XX
COMPLEX(EB) CIOR

!     .. ARRAY ARGUMENTS ..
REAL(EB)    RBIGA( * )
COMPLEX(EB) CBIGA( * )

!     .. LOCAL SCALARS ..
LOGICAL     DOWN
INTEGER     N
REAL(EB)    MIM, MRE, REZINV, RTMP
COMPLEX(EB) CTMP, ZINV

!                                  ** DECIDE WHETHER BIGA CAN BE
!                                  ** CALCULATED BY UP-RECURRENCE
MRE = REAL( CIOR )
MIM = ABS( AIMAG( CIOR ) )
IF( MRE<1.0_EB .OR. MRE>10.0_EB .OR. MIM>10.0_EB ) THEN
   DOWN = .TRUE.
ELSE IF( YESANG ) THEN
   DOWN = .TRUE.
!                                                    ** EQ. R48
   IF( MIM*XX < F2( MRE ) ) DOWN = .FALSE.
ELSE
   DOWN = .TRUE.
!                                                    ** EQ. R48
   IF( MIM*XX < F1( MRE ) ) DOWN = .FALSE.
END IF

ZINV   = 1.0_EB / ( CIOR*XX )
REZINV = 1.0_EB / ( MRE*XX )


IF( DOWN ) THEN
!                          ** COMPUTE INITIAL HIGH-ORDER BIGA USING
!                          ** LENTZ METHOD ( REF. 1, PP. 17-20 )
   CTMP = CONFRA( NTRM, ZINV )
!                                   *** DOWNWARD RECURRENCE FOR BIGA
   IF( NOABS ) THEN
!                                        ** NO-ABSORPTION CASE; EQ (R23)
      RBIGA( NTRM ) = REAL( CTMP )
      DO N = NTRM, 2, -1
         RBIGA( N - 1 ) = ( N*REZINV ) - 1.0_EB / ( ( N*REZINV ) + RBIGA( N ) )
      END DO
   ELSE
!                                         ** ABSORPTIVE CASE; EQ (R23)
      CBIGA( NTRM ) = CTMP
      DO N = NTRM, 2, -1
         CBIGA( N-1 ) = (N*ZINV) - 1.0_EB / ( (N*ZINV) + CBIGA( N ) )
   END DO
   END IF

ELSE
!                            *** UPWARD RECURRENCE FOR BIGA
   IF( NOABS ) THEN
!                                  ** NO-ABSORPTION CASE; EQ (R20,21)
      RTMP = SIN( MRE*XX )
      RBIGA( 1 ) = - REZINV + RTMP / ( RTMP*REZINV - COS( MRE*XX ) )

      DO N = 2, NTRM
         RBIGA( N ) = -( N*REZINV ) + 1.0_EB / ( ( N*REZINV ) - RBIGA( N - 1 ) )
      END DO

   ELSE
!                                     ** ABSORPTIVE CASE; EQ (R20,22)
      CTMP = EXP( - (0._EB,2._EB)*CIOR*XX )
      CBIGA( 1 ) = - ZINV + (1._EB-CTMP) /  ( ZINV * (1._EB-CTMP) - (0._EB,1._EB)*(1._EB+CTMP) )
      DO N = 2, NTRM
         CBIGA( N ) = - (N*ZINV) + 1.0_EB / ((N*ZINV) - CBIGA( N-1 ))
      END DO

END IF

END IF

CONTAINS

REAL(EB) FUNCTION F1(MRE)
REAL(EB), INTENT(IN) :: MRE
F1 = -8.0_EB + MRE**2*( 26.22_EB + MRE*( -0.4474_EB + MRE**3*( 0.00204_EB - 0.000175_EB*MRE ) ) )
END FUNCTION F1

REAL(EB) FUNCTION F2(MRE)  ! EQ. R47B
REAL(EB), INTENT(IN) :: MRE
F2 = 3.9_EB + MRE*( -10.8_EB + 13.78_EB*MRE )
END FUNCTION F2

END SUBROUTINE BIGA

COMPLEX(EB) FUNCTION CONFRA( N, ZINV )

!         COMPUTE BESSEL FUNCTION RATIO A-SUB-N FROM ITS
!         CONTINUED FRACTION USING LENTZ METHOD
!
!         ZINV = RECIPROCAL OF ARGUMENT OF A
!
!
!    I N T E R N A L    V A R I A B L E S
!    ------------------------------------
!
!    CAK      TERM IN CONTINUED FRACTION EXPANSION OF A (EQ. R25)
!    CAPT     FACTOR USED IN LENTZ ITERATION FOR A (EQ. R27)
!    CNUMER   NUMERATOR   IN CAPT  ( EQ. R28A )
!    CDENOM   DENOMINATOR IN CAPT  ( EQ. R28B )
!    CDTD     PRODUCT OF TWO SUCCESSIVE DENOMINATORS OF CAPT FACTORS
!                 ( EQ. R34C )
!    CNTN     PRODUCT OF TWO SUCCESSIVE NUMERATORS OF CAPT FACTORS
!                 ( EQ. R34B )
!    EPS1     ILL-CONDITIONING CRITERION
!    EPS2     CONVERGENCE CRITERION
!    KK       SUBSCRIPT K OF CAK  ( EQ. R25B )
!    KOUNT    ITERATION COUNTER ( USED TO PREVENT INFINITE LOOPING )
!    MAXIT    MAX. ALLOWED NO. OF ITERATIONS
!    MM       + 1  AND - 1, ALTERNATELY
! --------------------------------------------------------------------

!     .. SCALAR ARGUMENTS ..

INTEGER     N
COMPLEX(EB) ZINV

!     .. LOCAL SCALARS ..

INTEGER     KK, KOUNT, MAXIT, MM
REAL(EB)    EPS1, EPS2
COMPLEX(EB) CAK, CAPT, CDENOM, CDTD, CNTN, CNUMER
!
DATA      EPS1 / 1.E-2_EB / , EPS2 / 1.E-8_EB /
DATA      MAXIT / 10000 /
!                                 ** EQ. R25A
CONFRA = ( N + 1 ) * ZINV
MM     = - 1
KK     = 2*N + 3
!                                 ** EQ. R25B, K=2
CAK    = ( MM*KK ) * ZINV
CDENOM = CAK
CNUMER = CDENOM + 1.0_EB / CONFRA
KOUNT  = 1

L10: DO
   KOUNT = KOUNT + 1

   IF( KOUNT>MAXIT ) CALL ERRMSG('CONFRA--ITERATION FAILED TO CONVERGE',.TRUE.)

   MM  = - MM
   KK  = KK + 2
   !                                 ** EQ. R25B
   CAK = ( MM*KK ) * ZINV
   !                                          ** EQ. R32
   IF( ABS( CNUMER / CAK )<=EPS1 .OR. ABS( CDENOM / CAK )<=EPS1 ) THEN

   !                                  ** ILL-CONDITIONED CASE -- STRIDE
   !                                  ** TWO TERMS INSTEAD OF ONE

   !                                       ** EQ. R34
      CNTN   = CAK * CNUMER + 1.0_EB
      CDTD   = CAK * CDENOM + 1.0_EB
   !                                           ** EQ. R33
      CONFRA = ( CNTN / CDTD ) * CONFRA

      MM  = - MM
      KK  = KK + 2
   !                                 ** EQ. R25B
      CAK = ( MM*KK ) * ZINV
   !                                      ** EQ. R35
      CNUMER = CAK + CNUMER / CNTN
      CDENOM = CAK + CDENOM / CDTD
      KOUNT  = KOUNT + 1
      CYCLE L10

   ELSE
   !                           *** WELL-CONDITIONED CASE

   !                                  ** EQ. R27
      CAPT   = CNUMER / CDENOM
   !                                  ** EQ. R26
      CONFRA = CAPT * CONFRA
   !                                  ** CHECK FOR CONVERGENCE; EQ. R31

      IF (      ABS( REAL (CAPT) - 1.0_EB )>=EPS2 .OR. ABS( AIMAG(CAPT) )      >=EPS2 )  THEN
   !                                        ** EQ. R30
         CNUMER = CAK + 1.0_EB / CNUMER
         CDENOM = CAK + 1.0_EB / CDENOM
         CYCLE L10
      END IF
   END IF
   EXIT L10
END DO L10

RETURN
END FUNCTION CONFRA

SUBROUTINE MIPRNT( PRNT, XX, PERFCT, CREFIN, NUMANG, XMU, QEXT, QSCA, GQSC, NMOM, IPOLZN, MOMDIM, &
                   CALCMO, PMOM, SFORW, SBACK, TFORW, TBACK, S1, S2 )

!         PRINT SCATTERING QUANTITIES OF A SINGLE PARTICLE

!     .. SCALAR ARGUMENTS ..
LOGICAL     PERFCT
INTEGER     IPOLZN, MOMDIM, NMOM, NUMANG
REAL(EB)    GQSC, QEXT, QSCA, XX
COMPLEX(EB) CREFIN, SBACK, SFORW

!     .. ARRAY ARGUMENTS ..
LOGICAL     CALCMO( * ), PRNT( * )
REAL(EB)    PMOM( 0:MOMDIM, * ), XMU( * )
COMPLEX(EB) S1( * ), S2( * ), TBACK( * ), TFORW( * )

!     .. LOCAL SCALARS ..
CHARACTER   FMAT*22
INTEGER     I, J, M
REAL(EB)    FNORM, I1, I2

!     .. INTRINSIC FUNCTIONS ..
!     INTRINSIC AIMAG, CONJG, REAL


IF( PERFCT ) WRITE(LU_ERR, '(''1'',10X,A,1P,E11.4)' )  'PERFECTLY CONDUCTING CASE, SIZE PARAMETER =', XX

IF( .NOT.PERFCT ) WRITE(LU_ERR, '(''1'',10X,3(A,1P,E11.4))' ) &
    'REFRACTIVE INDEX:  REAL ', REAL( CREFIN ), '  IMAG ', &
    AIMAG( CREFIN ), ',   SIZE PARAMETER =', XX


IF( PRNT( 1 ) .AND. NUMANG>0 ) THEN

   WRITE(LU_ERR, '(/,A)' ) &
         '    COS(ANGLE)  ------- S1 ---------  ------- S2 ---------' &
         // '  --- S1*CONJG(S2) ---   I1=S1**2   I2=S2**2  (I1+I2)/2' &
         // '  DEG POLZN'

   DO I = 1, NUMANG
      I1 = REAL( S1( I ) )**2 + AIMAG( S1( I ) )**2
      I2 = REAL( S2( I ) )**2 + AIMAG( S2( I ) )**2
      WRITE(LU_ERR, '( I4, F10.6, 1P,10E11.3 )'   ) I, XMU(I), S1(I), S2(I), &
         S1(I)*CONJG(S2(I)),I1, I2, 0.5_EB*(I1+I2), (I2-I1)/(I2+I1)
   END DO
END IF


IF( PRNT( 2 ) ) THEN

   WRITE (LU_ERR, '(/,A,9X,A,17X,A,17X,A,/,(0P,F7.2, 1P,6E12.3) )' ) &
            '  ANGLE', 'S-SUB-1', 'T-SUB-1', 'T-SUB-2', &
                  0.0_EB,     SFORW,    TFORW(1),  TFORW(2), &
               180._EB,     SBACK,    TBACK(1),  TBACK(2)
   WRITE (LU_ERR, '(/,4(A,1P,E11.4))' ) &
            ' EFFICIENCY FACTORS,  EXTINCTION:', QEXT, &
                                 '   SCATTERING:', QSCA, &
                                 '   ABSORPTION:', QEXT-QSCA, &
                              '   RAD. PRESSURE:', QEXT-GQSC
   IF( NMOM>0 ) THEN
      WRITE(LU_ERR, '(/,A)' ) ' NORMALIZED MOMENTS OF :'
      IF( IPOLZN==0 ) WRITE(LU_ERR, '(''+'',27X,A)' ) 'PHASE FCN'
      IF( IPOLZN>0 ) WRITE(LU_ERR, '(''+'',33X,A)' )  'M1           M2          S21          D21'

      IF( IPOLZN<0 ) WRITE(LU_ERR, '(''+'',33X,A)' )  'R1           R2           R3           R4'
      FNORM = 4._EB/ ( XX**2 * QSCA )
      DO M = 0, NMOM
         WRITE(LU_ERR, '(A,I4)' ) '      MOMENT NO.', M
         DO J = 1, 4
            IF( CALCMO( J ) ) THEN
               WRITE( FMAT, '(A,I2,A)' ) '( ''+'', T', 24+(J-1)*13, ', 1P,E13.4 )'
               WRITE(LU_ERR, FMAT ) FNORM * PMOM( M, J )
            END IF
         END DO
      END DO
   END IF
END IF

RETURN
END SUBROUTINE MIPRNT

SUBROUTINE SMALL1( XX, NUMANG, XMU, QEXT, QSCA, GQSC, SFORW, SBACK, S1, S2, TFORW, TBACK, A, B )

!       SMALL-PARTICLE LIMIT OF MIE QUANTITIES IN TOTALLY REFLECTING
!       LIMIT ( MIE SERIES TRUNCATED AFTER 2 TERMS )
!
!        A,B       FIRST TWO MIE COEFFICIENTS, WITH NUMERATOR AND
!                  DENOMINATOR EXPANDED IN POWERS OF XX ( A FACTOR
!                  OF XX**3 IS MISSING BUT IS RESTORED BEFORE RETURN
!                  TO CALLING PROGRAM )  ( REF. 2, P. 1508 )

!     .. PARAMETERS ..
REAL(EB)  TWOTHR, FIVTHR, FIVNIN
PARAMETER ( TWOTHR = 2._EB/3._EB, FIVTHR = 5._EB/3._EB, FIVNIN = 5._EB/9._EB )

!     .. SCALAR ARGUMENTS ..
INTEGER     NUMANG
REAL(EB)    GQSC, QEXT, QSCA, XX
COMPLEX(EB) SBACK, SFORW

!     .. ARRAY ARGUMENTS ..
REAL(EB)    XMU( * )
COMPLEX(EB) A(*), B( * ), S1( * ), S2( * ), TBACK( * ), TFORW( * )

!     .. LOCAL SCALARS ..
INTEGER     J
REAL(EB)    RTMP

!                                                       ** EQ. R/A.5
A( 1 ) = CMPLX( 0._EB,    TWOTHR*( 1._EB - 0.2_EB*XX**2 ), EB ) / CMPLX( 1._EB - 0.5_EB*XX**2, TWOTHR*XX**3, EB )
!                                                      ** EQ. R/A.6
B( 1 ) = CMPLX( 0._EB, - ( 1._EB - 0.1_EB*XX**2 ) / 3._EB, EB ) / CMPLX( 1._EB + 0.5_EB*XX**2, - XX**3 / 3._EB, EB)
!                                                       ** EQ. R/A.7,8
A( 2 ) = CMPLX( 0._EB,   XX**2 / 30._EB, EB)
B( 2 ) = CMPLX( 0._EB, - XX**2 / 45._EB, EB)
!                                                       ** EQ. R/A.9
QSCA = 6._EB* XX**4 *( SQ( A(1) ) + SQ( B(1) ) + FIVTHR*( SQ( A(2) ) + SQ( B(2) ) ) )
QEXT = QSCA
!                                                       ** EQ. R/A.10
GQSC = 6._EB* XX**4 *REAL( A(1)*CONJG( A(2) + B(1) ) + ( B(1) + FIVNIN*A(2) )*CONJG( B(2) ) )

RTMP   = 1.5_EB * XX**3
SFORW  = RTMP*( A(1) + B(1) + FIVTHR*( A(2) + B(2) ) )
SBACK  = RTMP*( A(1) - B(1) - FIVTHR*( A(2) - B(2) ) )
TFORW( 1 ) = RTMP*( B(1) + FIVTHR*( 2._EB*B(2) - A(2) ) )
TFORW( 2 ) = RTMP*( A(1) + FIVTHR*( 2._EB*A(2) - B(2) ) )
TBACK( 1 ) = RTMP*( B(1) - FIVTHR*( 2._EB*B(2) + A(2) ) )
TBACK( 2 ) = RTMP*( A(1) - FIVTHR*( 2._EB*A(2) + B(2) ) )

DO J = 1, NUMANG
!                                                    ** EQ. R/A.11,12
   S1( J ) = RTMP*( A(1) + B(1)*XMU( J ) + FIVTHR*( A(2)*XMU( J ) +  B(2)*( 2._EB*XMU( J )**2 - 1._EB) ) )
   S2( J ) = RTMP*( B(1) + A(1)*XMU( J ) + FIVTHR*( B(2)*XMU( J ) +  A(2)*( 2._EB*XMU( J )**2 - 1._EB) ) )
END DO
!                                     ** RECOVER ACTUAL MIE COEFFICIENTS
A( 1 ) = XX**3 * A(1)
A( 2 ) = XX**3 * A(2)
B( 1 ) = XX**3 * B(1)
B( 2 ) = XX**3 * B(2)

RETURN
END SUBROUTINE SMALL1

SUBROUTINE SMALL2( XX, CIOR, CALCQE, NUMANG, XMU, QEXT, QSCA,GQSC, SFORW, SBACK, S1, S2, TFORW, TBACK,A, B )

!       SMALL-PARTICLE LIMIT OF MIE QUANTITIES FOR GENERAL REFRACTIVE
!       INDEX ( MIE SERIES TRUNCATED AFTER 2 TERMS )
!
!        A,B       FIRST TWO MIE COEFFICIENTS, WITH NUMERATOR AND
!                  DENOMINATOR EXPANDED IN POWERS OF XX ( A FACTOR
!                  OF XX**3 IS MISSING BUT IS RESTORED BEFORE RETURN
!                  TO CALLING PROGRAM )

!        CIORSQ    SQUARE OF REFRACTIVE INDEX

!     .. PARAMETERS ..
REAL(EB)    TWOTHR, FIVTHR
PARAMETER   ( TWOTHR = 2._EB/3._EB, FIVTHR = 5._EB/3._EB)

!     .. SCALAR ARGUMENTS ..
LOGICAL     CALCQE
INTEGER     NUMANG
REAL(EB)    GQSC, QEXT, QSCA, XX
COMPLEX(EB) CIOR, SBACK, SFORW

!     .. ARRAY ARGUMENTS ..
REAL(EB)    XMU( * )
COMPLEX(EB) A(*), B(*), S1(*), S2(*), TBACK(*), TFORW( * )

!     .. LOCAL SCALARS ..
INTEGER   J
REAL(EB)  RTMP
COMPLEX(EB)   CIORSQ, CTMP


CIORSQ = CIOR**2
CTMP   = CMPLX( 0._EB, TWOTHR, EB )*( CIORSQ - 1.0_EB )

!                                           ** EQ. R42A
A( 1 ) = CTMP*( 1._EB- 0.1_EB*XX**2 +   ( CIORSQ / 350._EB + 1._EB/280._EB)*XX**4 ) / &
          ( CIORSQ + 2._EB+ ( 1._EB- 0.7_EB*CIORSQ )*XX**2 -  &
          ( CIORSQ**2 / 175._EB- 0.275_EB*CIORSQ + 0.25_EB )*XX**4 + &
          XX**3 * CTMP * ( 1._EB- 0.1_EB*XX**2 ) )

!                                           ** EQ. R42B
B( 1 ) = ( XX**2 / 30._EB )*CTMP*( 1._EB+  ( CIORSQ / 35._EB - 1._EB/ 14._EB)*XX**2 ) /  &
         ( 1._EB- ( CIORSQ / 15._EB - 1._EB/6._EB)*XX**2 )

!                                           ** EQ. R42C

A( 2 ) = ( 0.1_EB*XX**2 )*CTMP*( 1._EB- XX**2 / 14._EB ) /  ( 2._EB*CIORSQ + 3._EB- &
         ( CIORSQ / 7._EB- 0.5_EB ) * XX**2 )

!                                           ** EQ. R40A

QSCA = (6._EB*XX**4) * ( SQ( A(1) ) + SQ( B(1) ) + FIVTHR * SQ( A(2) ) )

!                                           ** EQ. R40B
QEXT = QSCA
IF( CALCQE ) QEXT = 6._EB*XX * REAL( A(1) + B(1) + FIVTHR*A(2) )

!                                           ** EQ. R40C

GQSC = (6._EB*XX**4) * REAL( A(1)*CONJG( A(2) + B(1) ) )

RTMP   = 1.5_EB * XX**3
SFORW  = RTMP*( A(1) + B(1) + FIVTHR*A(2) )
SBACK  = RTMP*( A(1) - B(1) - FIVTHR*A(2) )
TFORW( 1 ) = RTMP*( B(1) - FIVTHR*A(2) )
TFORW( 2 ) = RTMP*( A(1) + 2._EB*FIVTHR*A(2) )
TBACK( 1 ) = TFORW(1)
TBACK( 2 ) = RTMP*( A(1) - 2._EB*FIVTHR*A(2) )


DO J = 1, NUMANG
!                                      ** EQ. R40D,E

   S1( J ) = RTMP*( A(1) + ( B(1) + FIVTHR*A(2) )*XMU( J ) )
   S2( J ) = RTMP*( B(1) + A(1)*XMU( J ) + FIVTHR*A(2)*( 2._EB*XMU( J )**2 - 1._EB) )
END DO

!                                     ** RECOVER ACTUAL MIE COEFFICIENTS
A( 1 ) = XX**3 * A(1)
A( 2 ) = XX**3 * A(2)
B( 1 ) = XX**3 * B(1)
B( 2 ) = ( 0._EB, 0._EB)

RETURN
END SUBROUTINE SMALL2

SUBROUTINE TESTMI( COMPAR, XX, CREFIN, MIMCUT, PERFCT, ANYANG, &
                   NMOM, IPOLZN, NUMANG, XMU, QEXT, QSCA, GQSC, &
                   SFORW, SBACK, S1, S2, TFORW, TBACK, PMOM, &
                   MOMDIM )

!         SET UP TO RUN TEST CASE WHEN  COMPAR = FALSE;  WHEN  = TRUE,
!         COMPARE MIE CODE TEST CASE RESULTS WITH CORRECT ANSWERS
!         AND ABORT IF EVEN ONE RESULT IS INACCURATE.

!         THE TEST CASE IS :  MIE SIZE PARAMETER = 10
!                             REFRACTIVE INDEX   = 1.5 - 0.1 I
!                             SCATTERING ANGLE = 140 DEGREES
!                             1 SEKERA MOMENT

!         RESULTS FOR THIS CASE MAY BE FOUND AMONG THE TEST CASES
!         AT THE END OF REFERENCE (1).

!         *** NOTE *** WHEN RUNNING ON SOME COMPUTERS, ESP. IN SINGLE
!         PRECISION, THE ACCUR CRITERION BELOW MAY HAVE TO BE RELAXED.
!         HOWEVER, IF ACCUR MUST BE SET LARGER THAN 10**-3 FOR SOME
!         SIZE PARAMETERS, YOUR COMPUTER IS PROBABLY NOT ACCURATE
!         ENOUGH TO DO MIE COMPUTATIONS FOR THOSE SIZE PARAMETERS.

!     ROUTINES CALLED :  ERRMSG, MIPRNT, TSTBAD

!     .. SCALAR ARGUMENTS ..

LOGICAL     ANYANG, COMPAR, PERFCT
INTEGER     IPOLZN, MOMDIM, NMOM, NUMANG
REAL(EB)    GQSC, MIMCUT, QEXT, QSCA, XX
COMPLEX(EB) CREFIN, SBACK, SFORW

!     .. ARRAY ARGUMENTS ..

REAL(EB)    PMOM( 0:MOMDIM, * ), XMU( * )
COMPLEX(EB) S1( * ), S2( * ), TBACK( * ), TFORW( * )

!     .. LOCAL SCALARS ..

LOGICAL     ANYSAV, OK, PERSAV
INTEGER     IPOSAV, M, N, NMOSAV, NUMSAV
REAL(EB)    MIMSAV, TESTGQ, TESTQE, TESTQS,XMUSAV, XXSAV
COMPLEX(EB) CRESAV, TESTS1, TESTS2, TESTSB, TESTSF

!     .. LOCAL ARRAYS ..

LOGICAL     CALCMO( 4 ), PRNT( 2 )
REAL(EB)    TESTPM( 0:1 )
COMPLEX(EB) TESTTB( 2 ), TESTTF( 2 )


SAVE      XXSAV, CRESAV, MIMSAV, PERSAV, ANYSAV, NMOSAV, IPOSAV,NUMSAV, XMUSAV

DATA      TESTQE / 2.459791_EB /,&
          TESTQS / 1.235144_EB /,&
          TESTGQ / 1.139235_EB /,&
          TESTSF / ( 61.49476_EB, -3.177994_EB ) /,&
          TESTSB / ( 1.493434_EB,  0.2963657_EB ) /,&
          TESTS1 / ( -0.1548380_EB, -1.128972_EB ) /,&
          TESTS2 / ( 0.05669755_EB, 0.5425681_EB ) /,&
          TESTTF / ( 12.95238_EB, -136.6436_EB ),&
                   ( 48.54238_EB, 133.4656_EB ) /,&
          TESTTB / ( 41.88414_EB, -15.57833_EB ),&
                   ( 43.37758_EB, -15.28196_EB ) /,&
          TESTPM / 227.1975_EB, 183.6898_EB /

IF( .NOT.COMPAR ) THEN
!                                   ** SAVE CERTAIN USER INPUT VALUES
   XXSAV  = XX
   CRESAV = CREFIN
   MIMSAV = MIMCUT
   PERSAV = PERFCT
   ANYSAV = ANYANG
   NMOSAV = NMOM
   IPOSAV = IPOLZN
   NUMSAV = NUMANG
   XMUSAV = XMU( 1 )
!                                   ** RESET INPUT VALUES FOR TEST CASE
   XX     = 10.0_EB
   CREFIN = ( 1.5_EB, -0.1_EB )
   MIMCUT = 0.0_EB
   PERFCT = .FALSE.
   ANYANG = .TRUE.
   NMOM   = 1
   IPOLZN = -1
   NUMANG = 1
   XMU( 1 ) = -0.7660444_EB

ELSE
!                                    ** COMPARE TEST CASE RESULTS WITH
!                                    ** CORRECT ANSWERS AND ABORT IF BAD
   OK = .TRUE.

   IF( WRONG( QEXT,TESTQE ) ) OK = TSTBAD( 'QEXT', ABS( ( QEXT - TESTQE ) / TESTQE ) )

   IF( WRONG( QSCA,TESTQS ) )OK = TSTBAD( 'QSCA', ABS( ( QSCA - TESTQS ) / TESTQS ) )

   IF( WRONG( GQSC,TESTGQ ) ) OK = TSTBAD( 'GQSC', ABS( ( GQSC - TESTGQ ) / TESTGQ ) )

   IF( WRONG( REAL( SFORW ),REAL( TESTSF ) ) .OR. WRONG( AIMAG( SFORW ),AIMAG( TESTSF ) ) ) &
         OK = TSTBAD( 'SFORW', ABS( ( SFORW - TESTSF ) / TESTSF ) )

   IF( WRONG( REAL( SBACK ),REAL( TESTSB ) ) .OR. WRONG( AIMAG( SBACK ),AIMAG( TESTSB ) ) ) &
         OK = TSTBAD( 'SBACK', ABS( ( SBACK - TESTSB ) / TESTSB ) )

   IF( WRONG( REAL( S1(1) ),REAL( TESTS1 ) ) .OR. WRONG( AIMAG( S1(1) ),AIMAG( TESTS1 ) ) ) &
         OK = TSTBAD( 'S1', ABS( ( S1(1) - TESTS1 ) / TESTS1 ) )

   IF( WRONG( REAL( S2(1) ),REAL( TESTS2 ) ) .OR. WRONG( AIMAG( S2(1) ),AIMAG( TESTS2 ) ) ) &
         OK = TSTBAD( 'S2', ABS( ( S2(1) - TESTS2 ) / TESTS2 ) )

   DO N = 1, 2
      IF( WRONG( REAL( TFORW(N) ),REAL( TESTTF(N) ) ) .OR. &
               WRONG( AIMAG( TFORW(N) ), &
               AIMAG( TESTTF(N) ) ) ) OK = TSTBAD( 'TFORW', &
               ABS( ( TFORW(N) - TESTTF(N) ) / TESTTF(N) ) )

      IF( WRONG( REAL( TBACK(N) ),REAL( TESTTB(N) ) ) .OR. &
               WRONG( AIMAG( TBACK(N) ), &
               AIMAG( TESTTB(N) ) ) ) OK = TSTBAD( 'TBACK', &
               ABS( ( TBACK(N) - TESTTB(N) ) / TESTTB(N) ) )
   END DO

   DO M = 0, 1
      IF ( WRONG( PMOM(M,1), TESTPM(M) ) ) OK =  TSTBAD( 'PMOM', ABS( (PMOM(M,1)-TESTPM(M)) / TESTPM(M) ) )
   END DO


   IF( .NOT.OK ) THEN
      PRNT( 1 ) = .TRUE.
      PRNT( 2 ) = .TRUE.
      CALCMO( 1 ) = .TRUE.
      CALCMO( 2 ) = .FALSE.
      CALCMO( 3 ) = .FALSE.
      CALCMO( 4 ) = .FALSE.
      CALL MIPRNT( PRNT, XX, PERFCT, CREFIN, NUMANG, XMU, QEXT, &
                        QSCA, GQSC, NMOM, IPOLZN, MOMDIM, CALCMO, PMOM, &
                        SFORW, SBACK, TFORW, TBACK, S1, S2 )
      CALL ERRMSG( 'MIEV0 -- SELF-TEST FAILED',.TRUE.)
   END IF
!                                       ** RESTORE USER INPUT VALUES
   XX     = XXSAV
   CREFIN = CRESAV
   MIMCUT = MIMSAV
   PERFCT = PERSAV
   ANYANG = ANYSAV
   NMOM   = NMOSAV
   IPOLZN = IPOSAV
   NUMANG = NUMSAV
   XMU( 1 ) = XMUSAV

END IF

CONTAINS

LOGICAL FUNCTION WRONG(CALC,EXACT)
REAL(EB) ACCUR
REAL(EB), INTENT(IN) :: CALC,EXACT
DATA  ACCUR / 1.E-4_EB /
WRONG = ABS( ( CALC - EXACT ) / EXACT )>ACCUR
END FUNCTION WRONG

END SUBROUTINE TESTMI


SUBROUTINE ERRMSG( MESSAG, FATAL )
! PRINT OUT A WARNING OR ERROR MESSAGE;  ABORT IF ERROR AFTER MAKING SYMBOLIC DUMP (MACHINE-SPECIFIC)
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
CHARACTER :: MESSAG*(*)
LOGICAL ::  FATAL,MSGLIM
INTEGER ::  MAXMSG, NUMMSG
SAVE      MAXMSG, NUMMSG, MSGLIM
DATA      NUMMSG / 0 /,  MAXMSG / 100 /,  MSGLIM /.FALSE./
IF (FATAL) CALL SHUTDOWN(MESSAG)
NUMMSG = NUMMSG + 1
IF( MSGLIM ) RETURN
IF( NUMMSG<=MAXMSG ) THEN
   WRITE(LU_ERR, '(/,2A,/)' ) ' ****** WARNING *****  ', MESSAG
ELSE
   WRITE(LU_ERR, 9000 )
   MSGLIM = .TRUE.
END IF
RETURN
9000 FORMAT( / , / , ' ****** TOO MANY WARNING MESSAGES --  ','THEY WILL NO LONGER BE PRINTED *******', / , / )
END SUBROUTINE ERRMSG


LOGICAL FUNCTION WRTBAD( VARNAM )
!          WRITE NAMES OF ERRONEOUS VARIABLES AND RETURN 'TRUE'
!      INPUT :   VARNAM = NAME OF ERRONEOUS VARIABLE TO BE WRITTEN
!                         ( CHARACTER, ANY LENGTH )
CHARACTER VARNAM*(*)
INTEGER   MAXMSG, NUMMSG
SAVE      NUMMSG, MAXMSG
DATA      NUMMSG / 0 /, MAXMSG / 50 /
WRTBAD = .TRUE.
NUMMSG = NUMMSG + 1
WRITE(LU_ERR, '(3A)' ) ' ****  INPUT VARIABLE  ', VARNAM,'  IN ERROR  ****'
IF( NUMMSG==MAXMSG ) CALL ERRMSG( 'TOO MANY INPUT ERRORS.  ABORTING...',.TRUE.)
RETURN
END FUNCTION WRTBAD


LOGICAL FUNCTION WRTDIM( DIMNAM, MINVAL )

! WRITE NAME OF TOO-SMALL SYMBOLIC DIMENSION AND THE VALUE IT SHOULD BE INCREASED TO;  RETURN 'TRUE'

!      INPUT :  DIMNAM = NAME OF SYMBOLIC DIMENSION WHICH IS TOO SMALL
!                        ( CHARACTER, ANY LENGTH )
!               MINVAL = VALUE TO WHICH THAT DIMENSION SHOULD BE
!                        INCREASED (AT LEAST)
CHARACTER :: DIMNAM*(*)
INTEGER ::  MINVAL

WRITE(LU_ERR, '(3A,I7)' ) ' ****  SYMBOLIC DIMENSION  ', DIMNAM,'  SHOULD BE INCREASED TO AT LEAST ', MINVAL
WRTDIM = .TRUE.
RETURN
END FUNCTION WRTDIM


LOGICAL FUNCTION TSTBAD( VARNAM, RELERR )
! WRITE NAME (VARNAM) OF VARIABLE FAILING SELF-TEST AND ITS PERCENT ERROR FROM THE CORRECT VALUE; RETURN 'FALSE'
CHARACTER :: VARNAM*(*)
REAL(EB) ::  RELERR
TSTBAD = .FALSE.
WRITE(LU_ERR, '(/,3A,1P,E11.2,A)' ) ' *** OUTPUT VARIABLE ', VARNAM, ' DIFFERED BY ', &
     100.*RELERR,' PER CENT FROM CORRECT VALUE.  SELF-TEST FAILED.'
RETURN
END FUNCTION TSTBAD


REAL(EB) FUNCTION SQ(CTMP)
COMPLEX(EB), INTENT(IN) :: CTMP
SQ = REAL( CTMP )**2 + AIMAG( CTMP )**2
END FUNCTION SQ

END MODULE MIEV

MODULE RAD

! Radiation heat transfer

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE OUTPUT_CLOCKS, ONLY: RADF_CLOCK,RADF_COUNTER
USE MESH_POINTERS
USE MESH_VARIABLES
USE RADCONS

IMPLICIT NONE (TYPE,EXTERNAL)
PRIVATE

PUBLIC INIT_RADIATION,COMPUTE_RADIATION,BLACKBODY_FRACTION

REAL(EB) :: TYY_FAC

CONTAINS


!> \brief Initialize radiation arrays.
!>
SUBROUTINE INIT_RADIATION

! Meanings of some variables defined here:
!
! THETAUP       Upper limit of spherical polar angle THETA on a solid angle interval
! THETALOW      Lower limit of spherical polar angle THETA on a solid angle interval
! PHIUP         Upper limit of spherical polar angle PHI on a solid angle interval
! PHILOW        Lower limit of spherical polar angle PHI on a solid angle interval
! PLANCK_C2     Second Planck radiation constant, 14387.69 micron.K (or 1.438769 cm.K)
! N, I, J, K    Integers used as indices of arrays in DO loops
! NRA           An integer that specifies the number of discrete radiation angles in the FVM
! NSB           Number of spectral bands in the piecewise constant absorption-emission spectrum
! IPC           An integer for looping over the Lagrangian particle classes


USE MEMORY_FUNCTIONS, ONLY : CHKMEMERR
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
USE MIEV
USE RADCAL_CALC
USE WSGG_ARRAYS
REAL(EB) :: THETAUP,THETALOW,PHIUP,PHILOW,F_THETA,PLANCK_C2,KSI,LT,RCRHO,YY,YY2,BBF,AP0,AMEAN,RADIANCE,TRANSMISSIVITY,X_N2,&
            THETA,PHI,DLO
INTEGER  :: N,I,J,K,IPC,IZERO,NN,NI,II,JJ,IIM,JJM,IBND,NS,NS2,NRA,NSB,RADCAL_TEMP(16)=0,RCT_SKIP=-1,IO
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
REAL(EB), ALLOCATABLE, DIMENSION(:) :: COSINE_ARRAY
TYPE (RAD_FILE_TYPE), POINTER :: RF

! A few miscellaneous constants

FOUR_SIGMA = 4._EB*SIGMA            ! Four times the Stefan-Boltzmann constant
RPI_SIGMA  = RPI*SIGMA              ! Stefan-Boltzmann constant divided by PI (RPI = reciprocal pi)

NRA = NUMBER_RADIATION_ANGLES
NSB = NUMBER_SPECTRAL_BANDS

! Set the opening angle of the cylindrical geometry equal to the azimuthal angle

IF (CYLINDRICAL) DPHI0 = PI/REAL(NRP(1))

ALLOCATE(RSA(1:NRA),STAT=IZERO)
CALL ChkMemErr('RADI','RSA',IZERO)
ALLOCATE(DLX(1:NRA),STAT=IZERO)
CALL ChkMemErr('RADI','DLX',IZERO)
ALLOCATE(DLY(1:NRA),STAT=IZERO)
CALL ChkMemErr('RADI','DLY',IZERO)
ALLOCATE(DLZ(1:NRA),STAT=IZERO)
CALL ChkMemErr('RADI','DLZ',IZERO)
IF (CYLINDRICAL) THEN
   ALLOCATE(DLB(1:NRA),STAT=IZERO)
   CALL ChkMemErr('RADI','DLB',IZERO)
ENDIF
ALLOCATE(DLN(-3:3,1:NRA),STAT=IZERO)
CALL ChkMemErr('RADI','DLN',IZERO)
ALLOCATE(DLM(1:NRA,3),STAT=IZERO)
CALL ChkMemErr('RADI','DLM',IZERO)
ALLOCATE(DLANG(3,1:NRA),STAT=IZERO)
CALL ChkMemErr('RADI','DLANG',IZERO)

! Determine mean direction normals and sweeping orders
! as described in the FDS Tech. Ref. Guide Vol. 1 Sec. 6.2.2.

N = 0
DO I=1,NRT
   DO J=1,NRP(I)
      N = N + 1
      THETALOW  = PI*REAL(I-1)/REAL(NRT)
      THETAUP   = PI*REAL(I)/REAL(NRT)
      F_THETA   = 0.5_EB*(THETAUP-THETALOW  - COS(THETAUP)*SIN(THETAUP) + COS(THETALOW)*SIN(THETALOW))
      THETA = 0.5_EB*(THETAUP+THETALOW)
      IF (CYLINDRICAL) THEN
         PHILOW = PI*REAL(J-1)/REAL(NRP(I))
         PHIUP  = PI*REAL(J)/REAL(NRP(I))
      ELSEIF (TWO_D) THEN
         PHILOW = TWOPI*REAL(J-1)/REAL(NRP(I)) + PIO2
         PHIUP  = TWOPI*REAL(J)/REAL(NRP(I))   + PIO2
      ELSE
         PHILOW = TWOPI*REAL(J-1)/REAL(NRP(I))
         PHIUP  = TWOPI*REAL(J)/REAL(NRP(I))
      ENDIF
      PHI=0.5_EB*(PHILOW+PHIUP)
      RSA(N) = (PHIUP-PHILOW)*(COS(THETALOW)-COS(THETAUP))
      IF (CYLINDRICAL) THEN
         DLX(N) =  (SIN(PHIUP)-SIN(PHILOW)) *F_THETA
         DLY(N) =  (-SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW))  +COS(DPHI0/2.)*(COS(PHILOW)-COS(PHIUP)))*F_THETA
         DLB(N) =  (-SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW))  -COS(DPHI0/2.)*(COS(PHILOW)-COS(PHIUP)))*F_THETA
         DLZ(N)    = 0.5_EB*(PHIUP-PHILOW)   * ((SIN(THETAUP))**2-(SIN(THETALOW))**2)
         DLANG(1,N)  = SIN(THETA)*COS(PHI)
         DLANG(2,N)  = SIN(THETA)*SIN(PHI)
         DLANG(3,N)  = COS(THETA)
         IF (N==1000000) WRITE(LU_ERR,'(A)') 'This line should never get executed. It is here only to prevent optimization.'
      ELSEIF (TWO_D) THEN
         DLX(N) = (SIN(PHIUP)-SIN(PHILOW))*F_THETA
         DLY(N) = 0._EB
         DLZ(N) = (COS(PHILOW)-COS(PHIUP))*F_THETA
         DLANG(1,N)  = COS(PHI)
         DLANG(2,N)  = 0._EB
         DLANG(3,N)  = SIN(PHI)
      ELSE
         DLX(N) = (SIN(PHIUP)-SIN(PHILOW))*F_THETA
         DLY(N) = (COS(PHILOW)-COS(PHIUP))*F_THETA
         DLZ(N)    = 0.5_EB*(PHIUP-PHILOW)      * ((SIN(THETAUP))**2-(SIN(THETALOW))**2)
         DLANG(1,N)  = SIN(THETA)*COS(PHI)
         DLANG(2,N)  = SIN(THETA)*SIN(PHI)
         DLANG(3,N)  = COS(THETA)
      ENDIF
   ENDDO
ENDDO

! Set (wall normal)*(angle vector) value

DO N = 1,NRA
   DLN( 0,N) = 0._EB !prevent undefined variable errors
   DLN(-1,N) = -DLX(N)
   DLN( 1,N) =  DLX(N)
   DLN(-2,N) = -DLY(N)
   DLN( 2,N) =  DLY(N)
   DLN(-3,N) = -DLZ(N)
   DLN( 3,N) =  DLZ(N)
ENDDO

! In axially symmetric case, each angle represents two symmetric angles. So weight the intensities by two.

WEIGH_CYL = 1._EB
IF (CYLINDRICAL) THEN
   WEIGH_CYL = 2._EB
   ! Wall direction cosines are only used for flux integrations, so they can by multiplied in advance.
   DLN = WEIGH_CYL * DLN
ENDIF

! Calculate mirroring matrix

N = 0
DO I=1,NRT
   DO J=1,NRP(I)
      N = N + 1
      DO K=1,3
         IF (TWO_D .AND. .NOT.CYLINDRICAL) THEN
            SELECT CASE(K)
               CASE(1)             ! X-surfaces
                  IIM = 1
                  JJM = NRP(I) - J + 1
               CASE(2)             ! Y-surfaces
                  IIM = 1
                  JJM = J
               CASE(3)             ! Z-surfaces
                  IIM = 1
                  JJM = NRP(I)/2 - J + 1
            END SELECT
            JJM = MODULO(JJM,NRP(I))
            IF (JJM==0) JJM = NRP(I)
         ELSE
            SELECT CASE(K)
               CASE(1)             ! X-surfaces
                  IIM = I
                  JJM = NRP(I)/2 - J + 1
               CASE(2)             ! Y-surfaces
                  IIM = I
                  JJM = NRP(I) - J + 1
               CASE(3)             ! Z-surfaces
                  IIM = NRT - I + 1
                  JJM = J
            END SELECT
            IIM = MODULO(IIM,NRT)
            JJM = MODULO(JJM,NRP(I))
            IF (IIM==0) IIM = NRT
            IF (JJM==0) JJM = NRP(I)
         ENDIF

         NN = 0
         DO II = 1,IIM
            DO JJ = 1,NRP(II)
               NN = NN + 1
               IF ((II==IIM).AND.(JJ==JJM)) NI = NN
            ENDDO
         ENDDO
         DLM(N,K) = NI
      ENDDO
   ENDDO
ENDDO

!-----------------------------------------------------
!
!            Radiative properties computation
!
!-----------------------------------------------------


! General parameters

RTMPMAX = 2470._EB     ! Maximum temperature for property tables
RTMPMIN = 270._EB      ! Minimum temperature for property tables

! Setup spectral information

INIT_WIDE_BAND: IF (WIDE_BAND_MODEL) THEN

   ! Fraction of blackbody emission in a wavelength interval

   PLANCK_C2 = 14387.69_EB       ! Value of the 2nd Planck radiation constant in micron.K
   NLAMBDAT  = 4000
   LTSTEP    = 25.0_EB           ! maximum LAMBDA*T = NLANBDAT*LTSTEP
   ALLOCATE(BBFRAC(0:NLAMBDAT),STAT=IZERO)
   CALL ChkMemErr('INIT','BBFRAC',IZERO)

   BBFRAC = 0._EB
   LT     = 0._EB

   DO I = 1,NLAMBDAT
      LT  = LT + LTSTEP
      KSI = PLANCK_C2/LT

      DO J = 1,50
         BBFRAC(I) = BBFRAC(I) + EXP(-KSI*REAL(J))/REAL(J) * (KSI**3 + 3.*KSI**2/REAL(J) + 6.*KSI/REAL(J)**2 + 6./REAL(J)**3)
      ENDDO
   ENDDO

   BBFRAC =  BBFRAC * 15._EB/PI**4

   ! Define band limit wave lengths in micrometers

   IF (.NOT.ALLOCATED(WL_LOW).OR. .NOT.ALLOCATED(WL_HIGH)) THEN

      ALLOCATE(WL_LOW(1:NSB),STAT=IZERO)
      CALL ChkMemErr('INIT','WL_LOW',IZERO)
      ALLOCATE(WL_HIGH(1:NSB),STAT=IZERO)
      CALL ChkMemErr('INIT','WL_HIGH',IZERO)

      ! Define the band limits as function of the fuel used
      ! Use (SPECIES(FUEL_INDEX)%RADCAL_ID (FUEL_INDEX global variable defined in module GLOBAL_CONSTANTS

      IF (FUEL_INDEX>0) THEN

         SELECT CASE(SPECIES(FUEL_INDEX)%RADCAL_ID)
            CASE('METHANE')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.630_EB, 2.940_EB, 4.170_EB, 4.600_EB,  10.000_EB /)
               WL_HIGH(1:NSB) = (/ 2.630_EB, 2.940_EB, 4.170_EB, 4.600_EB, 10.00_EB, 200.000_EB /)
            CASE('ETHANE')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.632_EB, 2.985_EB, 4.000_EB, 6.061_EB,   9.174_EB /)
               WL_HIGH(1:NSB) = (/ 2.632_EB, 2.985_EB, 4.000_EB, 6.061_EB, 9.174_EB, 200.000_EB /)
            CASE('ETHYLENE')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.632_EB, 2.963_EB, 3.571_EB, 6.061_EB,  12.821_EB /)
               WL_HIGH(1:NSB) = (/ 2.632_EB, 2.963_EB, 3.571_EB, 6.061_EB,12.821_EB, 200.000_EB /)
            CASE('PROPANE')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.632_EB, 2.985_EB, 3.922_EB, 6.061_EB,   8.511_EB /)
               WL_HIGH(1:NSB) = (/ 2.632_EB, 2.985_EB, 3.922_EB, 6.061_EB, 8.511_EB, 200.000_EB /)
            CASE('PROPYLENE')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.632_EB, 3.077_EB, 3.846_EB, 5.128_EB,   8.511_EB /)
               WL_HIGH(1:NSB) = (/ 2.632_EB, 3.077_EB, 3.846_EB, 5.128_EB, 8.511_EB, 200.000_EB /)
            CASE('N-HEPTANE')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.632_EB, 3.077_EB, 3.922_EB, 5.634_EB,   9.091_EB /)
               WL_HIGH(1:NSB) = (/ 2.632_EB, 3.077_EB, 3.922_EB, 5.634_EB, 9.091_EB, 200.000_EB /)
            CASE('TOLUENE')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.632_EB, 3.125_EB, 3.922_EB, 4.878_EB,   8.333_EB /)
               WL_HIGH(1:NSB) = (/ 2.632_EB, 3.125_EB, 3.922_EB, 4.878_EB, 8.333_EB, 200.000_EB /)
            CASE('METHANOL')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.614_EB, 3.125_EB, 3.846_EB, 5.970_EB,   8.889_EB /)
               WL_HIGH(1:NSB) = (/ 2.614_EB, 3.125_EB, 3.846_EB, 5.970_EB, 8.889_EB, 200.000_EB /)
            CASE('MMA')
               WL_LOW(1:NSB)  = (/ 1.000_EB, 2.632_EB, 3.077_EB, 3.774_EB, 4.878_EB,   8.000_EB /)
               WL_HIGH(1:NSB) = (/ 2.632_EB, 3.077_EB, 3.774_EB, 4.878_EB, 8.000_EB, 200.000_EB /)
         END SELECT

      ELSE

         ! Use the methane bands if there is no fuel species

         WL_LOW(1:NSB)  = (/ 1.000_EB, 2.630_EB, 2.940_EB, 4.170_EB, 4.600_EB,  10.000_EB /)
         WL_HIGH(1:NSB) = (/ 2.630_EB, 2.940_EB, 4.170_EB, 4.600_EB, 10.00_EB, 200.000_EB /)

      ENDIF

   ENDIF

ENDIF INIT_WIDE_BAND

!----------------------------------------------------------------------------
!
!     Tables for gas phase absorption coefficient
!
!     CONTROLLING PROGRAM FOR SUBROUTINE "RADCAL", A NARROW-BAND
!     MODEL FOR CALCULATING SPECTRAL INTENSITY (W/M-2/SR/MICRON) AND
!     SPECTRAL TRANSMITTANCE VERSUS WAVELENGTH (MICRONS) IN A NONISO-
!     THERMAL, VARIABLE COMPOSITION  MIXTURE OF CO2, H2O, CO, N2, O2,
!     CH4, AND SOOT. FOR A HOMOGENEOUS PATH, THE PROGRAM ALSO COMPUTES
!     THE PLANCK-MEAN ABSORPTION COEF, AP0, THE INCIDENT-MEAN ABSORPTION
!     COEFFICIENT, AIWALL, AND THE EFFECTIVE-MEAN ABSORPTION COEFFICIENT,
!     AMEAN, ALL IN UNITS OF INVERSE METERS.
!
!     INPUT PARAMETERS:
!          NPT=NUMBER OF HOMOGENEOUS ELEMENTS
!          SEGMENT_LENGTH_M(J)=THICKNESS OF J TH ELEMENT, M
!          TEMP_GAS(J)=TEMPERATURE OF J TH ELEMENT, K.
!          PARTIAL_PRESSURES_ATM(I,J)=PARTIAL PRESSURE OF GASEOUS COMPONENTS I, kPa:
!          OMMIN=MINIMUM WAVE NUMBER IN SPECTRUM, CM-1.
!          OMMAX=MAXIMUM WAVE NUMBER IN SPECTRUM, CM-1.
!
!-------------------------------------------------------------------------

MAKE_KAPPA_ARRAYS: IF (.NOT.SOLID_PHASE_ONLY .AND. ANY(SPECIES%RADCAL_ID/='null') .AND. .NOT.WSGG_MODEL) THEN

   ! Check for valid RADCAL species and setup arrays from ZZ to RADCAL_YY

   N_RADCAL_ARRAY_SIZE = 0
   GET_RADCAL_SPECIES: DO NS=1,N_SPECIES
      SELECT CASE (SPECIES(NS)%RADCAL_ID)
         CASE('CARBON DIOXIDE')
            IF (RADCAL_TEMP(1)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=1
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='CARBON DIOXIDE'
               RADCAL_TEMP(1)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(1)
         CASE('WATER VAPOR')
            IF (RADCAL_TEMP(2)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=2
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='WATER VAPOR'
               SPECIES(NS)%RADCAL_INDEX=N_RADCAL_ARRAY_SIZE
               RADCAL_TEMP(2)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(2)
         CASE('CARBON MONOXIDE')
            IF (RADCAL_TEMP(3)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=3
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='CARBON MONOXIDE'
               RADCAL_TEMP(3)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(3)
         CASE('METHANE')
            IF (RADCAL_TEMP(4)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=4
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='METHANE'
               RADCAL_TEMP(4)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(4)
         CASE('ETHYLENE')
            IF (RADCAL_TEMP(5)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=5
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='ETHYLENE'
               RADCAL_TEMP(5)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(5)
         CASE('ETHANE')
            IF (RADCAL_TEMP(6)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=6
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='ETHANE'
               RADCAL_TEMP(6)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(6)
         CASE('PROPYLENE')
            IF (RADCAL_TEMP(7)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=7
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='PROPYLENE'
               RADCAL_TEMP(7)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(7)
         CASE('PROPANE')
            IF (RADCAL_TEMP(8)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=8
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='PROPANE'
               RADCAL_TEMP(8)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(8)
         CASE('TOLUENE')
            IF (RADCAL_TEMP(9)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=9
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='TOLUENE'
               RADCAL_TEMP(9)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(9)
         CASE('N-HEPTANE')
            IF (RADCAL_TEMP(10)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=10
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='N-HEPTANE'
               RADCAL_TEMP(10)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(10)
         CASE('METHANOL')
            IF (RADCAL_TEMP(11)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=11
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='METHANOL'
               RADCAL_TEMP(11)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(11)
         CASE('MMA')
            IF (RADCAL_TEMP(12)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=12
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='MMA'
               RADCAL_TEMP(12)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(12)
         CASE('SOOT')
            IF (RADCAL_TEMP(16)==0) THEN
               N_RADCAL_ARRAY_SIZE=N_RADCAL_ARRAY_SIZE+1
               RADCAL_SPECIES_INDEX(N_RADCAL_ARRAY_SIZE)=16
               RADCAL_SPECIES_ID(N_RADCAL_ARRAY_SIZE)='SOOT'
               RADCAL_TEMP(16)=N_RADCAL_ARRAY_SIZE
            ENDIF
            SPECIES(NS)%RADCAL_INDEX=RADCAL_TEMP(16)
      END SELECT
   END DO GET_RADCAL_SPECIES

   BUILD_KAPPA_ARRAY: IF (N_RADCAL_ARRAY_SIZE>0) THEN

      ALLOCATE(Z2RADCAL_SPECIES(N_RADCAL_ARRAY_SIZE,1:N_TRACKED_SPECIES),STAT=IZERO)
      CALL ChkMemErr('RADI','ZZ2RADCAL_SPECIES',IZERO)
      Z2RADCAL_SPECIES = 0._EB

      DO NS=1,N_TRACKED_SPECIES
         IF (SPECIES_MIXTURE(NS)%EVAPORATION_SMIX_INDEX>0) THEN
            CALL MEAN_CROSS_SECTIONS(SMIX_INDEX=NS)
            CYCLE !No gas absorption for the liquid phase
         ENDIF
         DO NS2=1,N_SPECIES
            IF (SPECIES(NS2)%RADCAL_INDEX > 0) THEN
               IF (SPECIES(NS2)%RADCAL_ID/='SOOT') THEN
                  Z2RADCAL_SPECIES(SPECIES(NS2)%RADCAL_INDEX,NS) = Z2RADCAL_SPECIES(SPECIES(NS2)%RADCAL_INDEX,NS) + &
                                                                   REAL(N_KAPPA_Y,EB)**4 / SPECIES(NS2)%MW * Z2Y(NS2,NS)
               ELSE
                  Z2RADCAL_SPECIES(SPECIES(NS2)%RADCAL_INDEX,NS) = Z2RADCAL_SPECIES(SPECIES(NS2)%RADCAL_INDEX,NS) + &
                                                                   REAL(N_KAPPA_Y,EB)**4 * 5._EB * Z2Y(NS2,NS) * &
                                                                   SPECIES(SOOT_INDEX)%DENSITY_SOLID/SPECIES(NS2)%DENSITY_SOLID
               ENDIF
            ENDIF
         ENDDO
      ENDDO

      ! Allocate arrays for RadCal

      CALL RCALLOC

      ! Set the Mean Beam Length to 10 cm unless the user desires otherwise

      IF (PATH_LENGTH < 0._EB) PATH_LENGTH = 0.1_EB
      ALLOCATE(SEGMENT_LENGTH_M(1))
      ALLOCATE(TOTAL_PRESSURE_ATM(1))
      ALLOCATE(TEMP_GAS(1))
      ALLOCATE(PARTIAL_PRESSURES_ATM(16,1))
      SEGMENT_LENGTH_M(1) = MAX(PATH_LENGTH,1.0E-4_EB)
      TOTAL_PRESSURE_ATM(1)=P_INF/P_STP
      NPT = 1
      TWALL = RADTMP
      LAMBDAMIN = -1.1E+4_EB
      LAMBDAMAX = -1.0E+4_EB

      ! Using RadCal, create look-up tables for the absorption coefficients for all gas species, mixture fraction or aerosols

      ALLOCATE (RADCAL_SPECIES2KAPPA(N_RADCAL_ARRAY_SIZE,0:N_KAPPA_Y,0:N_KAPPA_T,NSB),STAT=IZERO)
      CALL ChkMemErr('RADI','RADCAL_SPECIES2KAPPA',IZERO)
      RADCAL_SPECIES2KAPPA = 0._EB
      BBF = 1._EB
      OMMIN = 50._EB
      OMMAX = 10000._EB
      BAND_LOOP_Z: DO IBND = 1,NSB
         IF (NSB>1) THEN
            OMMIN = REAL(NINT(1.E4_EB/WL_HIGH(IBND)),EB)
            OMMAX = REAL(NINT(1.E4_EB/WL_LOW(IBND)),EB)
         ENDIF
         CALL INIT_RADCAL
         T_LOOP_Z: DO K = 0,N_KAPPA_T
            TEMP_GAS(1) = RTMPMIN + K*(RTMPMAX-RTMPMIN)/N_KAPPA_T
            ! AMEAN will not be calculated close to RADTMP, where it cannot be solved
            IF (ABS(TEMP_GAS(1)-RADTMP)<=0.4_EB*(RTMPMAX-RTMPMIN)/N_KAPPA_T .AND. PATH_LENGTH > 0.0_EB) THEN
               RCT_SKIP = K
               CYCLE T_LOOP_Z
            ENDIF
            RCRHO = MW_AIR*P_INF/(R0*TEMP_GAS(1))
            IF (NSB>1) BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TEMP_GAS(1))
            Y_LOOP_Z: DO J=1,N_KAPPA_Y
               YY = (REAL(J,EB)/REAL(N_KAPPA_Y,EB))**4
               YY2 = 1._EB-YY
               X_N2 = YY2/28._EB
               N = 0
               RADCAL_SPECIES_LOOP: DO NS = 1, N_RADCAL_ARRAY_SIZE
                  PARTIAL_PRESSURES_ATM = 0._EB
                  SELECT CASE(RADCAL_SPECIES_INDEX(NS))
                     CASE(1) ! CARBON DIOXIDE
                        PARTIAL_PRESSURES_ATM(1,1) = YY/(YY+44._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(1,1)
                     CASE(2) ! WATER VAPOR
                        PARTIAL_PRESSURES_ATM(2,1) = YY/(YY+18._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(2,1)
                     CASE(3) ! CARBON MONOXIDE
                        PARTIAL_PRESSURES_ATM(3,1) = YY/(YY+28._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(3,1)
                     CASE(4) ! METHANE
                        PARTIAL_PRESSURES_ATM(4,1) = YY/(YY+16._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(4,1)
                     CASE(5) ! EHTYLENE
                        PARTIAL_PRESSURES_ATM(5,1) = YY/(YY+28._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(5,1)
                     CASE(6) ! ETHANE
                        PARTIAL_PRESSURES_ATM(6,1) = YY/(YY+30._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(6,1)
                     CASE(7) ! PROPYLENE
                        PARTIAL_PRESSURES_ATM(7,1) = YY/(YY+42._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(7,1)
                     CASE(8) ! PROPANE
                        PARTIAL_PRESSURES_ATM(8,1) = YY/(YY+44._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(8,1)
                     CASE(9) ! TOLUENE
                        PARTIAL_PRESSURES_ATM(9,1) = YY/(YY+92._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(9,1)
                     CASE(10) ! N-HEPTANE
                        PARTIAL_PRESSURES_ATM(10,1) = YY/(YY+100._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(10,1)
                     CASE(11) ! METHANOL
                        PARTIAL_PRESSURES_ATM(11,1) = YY/(YY+32._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(11,1)
                     CASE(12) ! MMA
                        PARTIAL_PRESSURES_ATM(12,1) = YY/(YY+100._EB*X_N2)
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB - PARTIAL_PRESSURES_ATM(12,1)
                     CASE(16) ! SOOT
                        YY2 = 0.2_EB*YY
                        PARTIAL_PRESSURES_ATM(16,1) = YY2*RCRHO/SPECIES(SOOT_INDEX)%DENSITY_SOLID
                        PARTIAL_PRESSURES_ATM(14,1) = 1._EB
                  END SELECT
                  CALL SUB_RADCAL(AMEAN,AP0,RADIANCE,TRANSMISSIVITY)
                  IF (NSB==1 .AND. PATH_LENGTH > 0.0_EB) THEN
                     RADCAL_SPECIES2KAPPA(NS,J,K,1) = MIN(AMEAN,AP0)
                  ELSE
                     RADCAL_SPECIES2KAPPA(NS,J,K,IBND) = AMEAN/BBF
                  ENDIF
               END DO RADCAL_SPECIES_LOOP
            ENDDO Y_LOOP_Z
         ENDDO T_LOOP_Z
         ! Interpolate KAPPA at RADTMP
         IF (RCT_SKIP == 0) THEN
            RADCAL_SPECIES2KAPPA(:,:,0,IBND) = RADCAL_SPECIES2KAPPA(:,:,1,IBND)
         ELSEIF (RCT_SKIP == N_KAPPA_T) THEN
            RADCAL_SPECIES2KAPPA(:,:,N_KAPPA_T,IBND) = RADCAL_SPECIES2KAPPA(:,:,N_KAPPA_T-1,IBND)
         ELSEIF (RCT_SKIP > 0) THEN
            RADCAL_SPECIES2KAPPA(:,:,RCT_SKIP,IBND) = 0.5_EB*(RADCAL_SPECIES2KAPPA(:,:,RCT_SKIP-1,IBND)+ &
                                                              RADCAL_SPECIES2KAPPA(:,:,RCT_SKIP+1,IBND))
         ENDIF
         CALL RCDEALLOC2  ! Deallocate RadCal wavelength dependent arrays
      ENDDO BAND_LOOP_Z
      !Adjust values from /cm to /m
      RADCAL_SPECIES2KAPPA =  RADCAL_SPECIES2KAPPA * 100._EB
      CALL RCDEALLOC  ! Deallocate RadCal arrays

   ENDIF BUILD_KAPPA_ARRAY

   ! Trap any errors

   IF (ANY(RADCAL_SPECIES2KAPPA<0._EB)) CALL SHUTDOWN('ERROR: KAPPA < 0 in RADCAL')

ENDIF MAKE_KAPPA_ARRAYS

! Reference: Bordbar, Wecel, Hyppanen.  A line by line based weighted sum of gray gases model for inhomogeneous
!            CO2-H2O mixture in oxy-fired comustion.  Combustion and Flame, 161 (2014) 2435-2445.

MAKE_WSGG_ARRAYS: IF (.NOT.SOLID_PHASE_ONLY .AND. WSGG_MODEL) THEN
   ALLOCATE(WSGG_B1_ARRAY(1:4,0:4))
   ALLOCATE(WSGG_B2_ARRAY(1:4,0:4))
   ALLOCATE(WSGG_C_ARRAY(1:4,0:4,0:4))
   ALLOCATE(WSGG_D_ARRAY(1:4,0:4))
   ALLOCATE(WSGG_KAPPAP1_ARRAY(1:5))
   ALLOCATE(WSGG_KAPPAP2_ARRAY(1:5))

   WSGG_B1_ARRAY(1,0:4) = (/  8.425766E-1_EB, -1.442229E+0_EB,  1.286974E+0_EB, -5.202712E-1_EB,  7.581559E-2_EB /)
   WSGG_B1_ARRAY(2,0:4) = (/ -3.023864E-2_EB,  5.264245E-1_EB, -6.209696E-1_EB,  2.704755E-1_EB, -4.090690E-2_EB /)
   WSGG_B1_ARRAY(3,0:4) = (/  1.070243E-1_EB, -1.989596E-1_EB,  3.101602E-1_EB, -1.737230E-1_EB,  3.081180E-2_EB /)
   WSGG_B1_ARRAY(4,0:4) = (/  3.108972E-2_EB,  1.981489E-1_EB, -2.543676E-1_EB,  1.061331E-1_EB, -1.498231E-2_EB /)

   WSGG_B2_ARRAY(1,0:4) = (/  7.129509E-1_EB, -1.378353E+0_EB,  1.555028E+0_EB, -6.636291E-1_EB,  9.773674E-2_EB /)
   WSGG_B2_ARRAY(2,0:4) = (/  1.589917E-1_EB,  5.635578E-2_EB,  2.666874E-1_EB, -2.040335E-1_EB,  3.742408E-2_EB /)
   WSGG_B2_ARRAY(3,0:4) = (/ -1.196373E-1_EB,  1.349665E+0_EB, -1.544797E+0_EB,  6.397595E-1_EB, -9.153650E-2_EB /)
   WSGG_B2_ARRAY(4,0:4) = (/  3.078250E-1_EB, -6.003555E-1_EB,  4.441261E-1_EB, -1.468813E-1_EB,  1.824702E-2_EB /)

   WSGG_C_ARRAY(1,0,0:4) = (/  0.7412956_EB, -0.5244441_EB,  0.5822860_EB, -0.2096994_EB,  0.0242031_EB /)
   WSGG_C_ARRAY(1,1,0:4) = (/ -0.9412652_EB,  0.2799577_EB, -0.7672319_EB,  0.3204027_EB, -0.0391017_EB /)
   WSGG_C_ARRAY(1,2,0:4) = (/  0.8531866_EB,  0.0823075_EB,  0.5289430_EB, -0.2468463_EB,  0.0310940_EB /)
   WSGG_C_ARRAY(1,3,0:4) = (/ -0.3342806_EB,  0.1474987_EB, -0.4160689_EB,  0.1697627_EB, -0.0204066_EB /)
   WSGG_C_ARRAY(1,4,0:4) = (/  0.0431436_EB, -0.0688622_EB,  0.1109773_EB, -0.0420861_EB,  0.0049188_EB /)
   WSGG_C_ARRAY(2,0,0:4) = (/  0.1552073_EB, -0.4862117_EB,  0.3668088_EB, -0.1055508_EB,  0.0105857_EB /)
   WSGG_C_ARRAY(2,1,0:4) = (/  0.6755648_EB,  1.4092710_EB, -1.3834490_EB,  0.4575210_EB, -0.0501976_EB /)
   WSGG_C_ARRAY(2,2,0:4) = (/ -1.1253940_EB, -0.5913199_EB,  0.9085441_EB, -0.3334201_EB,  0.0384236_EB /)
   WSGG_C_ARRAY(2,3,0:4) = (/  0.6040543_EB, -0.0553385_EB, -0.1733014_EB,  0.0791608_EB, -0.0098934_EB /)
   WSGG_C_ARRAY(2,4,0:4) = (/ -0.1105453_EB,  0.0464663_EB, -0.0016129_EB, -0.0035398_EB,  0.0006121_EB /)
   WSGG_C_ARRAY(3,0,0:4) = (/  0.2550242_EB,  0.3805403_EB, -0.4249709_EB,  0.1429446_EB, -0.0157408_EB /)
   WSGG_C_ARRAY(3,1,0:4) = (/ -0.6065428_EB,  0.3494024_EB,  0.1853509_EB, -0.1013694_EB,  0.0130244_EB /)
   WSGG_C_ARRAY(3,2,0:4) = (/  0.8123855_EB, -1.1020090_EB,  0.4046178_EB, -0.0811822_EB,  0.0062981_EB /)
   WSGG_C_ARRAY(3,3,0:4) = (/ -0.4532290_EB,  0.6784475_EB, -0.3432603_EB,  0.0883088_EB, -0.0084152_EB /)
   WSGG_C_ARRAY(3,4,0:4) = (/  0.0869309_EB, -0.1306996_EB,  0.0741446_EB, -0.0202929_EB,  0.0020110_EB /)
   WSGG_C_ARRAY(4,0,0:4) = (/ -0.0345199_EB,  0.2656726_EB, -0.1225365_EB,  0.0300151_EB, -0.0028205_EB /)
   WSGG_C_ARRAY(4,1,0:4) = (/  0.4112046_EB, -0.5728350_EB,  0.2924490_EB, -0.0798076_EB,  0.0079966_EB /)
   WSGG_C_ARRAY(4,2,0:4) = (/ -0.5055995_EB,  0.4579559_EB, -0.2616436_EB,  0.0764841_EB, -0.0079084_EB /)
   WSGG_C_ARRAY(4,3,0:4) = (/  0.2317509_EB, -0.1656759_EB,  0.1052608_EB, -0.0321935_EB,  0.0033870_EB /)
   WSGG_C_ARRAY(4,4,0:4) = (/ -0.0375491_EB,  0.0229520_EB, -0.0160047_EB,  0.0050463_EB, -0.0005364_EB /)

   WSGG_KAPPAP1_ARRAY(1:5) = (/ 3.388079E-2_EB, 4.544269E-1_EB, 4.680226_EB, 1.038439E2_EB, 0._EB /)

   WSGG_KAPPAP2_ARRAY(1:5) = (/ 7.703541E-2_EB, 8.242941E-1_EB, 6.854761_EB, 6.593653E1_EB, 0._EB /)

   WSGG_D_ARRAY(1,0:4) = (/ 0.0340429_EB,  0.0652305_EB, -0.0463685_EB,  0.0138684_EB, -0.0014450_EB /)
   WSGG_D_ARRAY(2,0:4) = (/ 0.3509457_EB,  0.7465138_EB, -0.5293090_EB,  0.1594423_EB, -0.0166326_EB /)
   WSGG_D_ARRAY(3,0:4) = (/ 4.5707400_EB,  2.1680670_EB, -1.4989010_EB,  0.4917165_EB, -0.0542999_EB /)
   WSGG_D_ARRAY(4,0:4) = (/ 109.81690_EB, -50.923590_EB,  23.432360_EB, -5.1638920_EB,  0.4393889_EB /)
ENDIF MAKE_WSGG_ARRAYS

TYY_FAC=N_KAPPA_T / (RTMPMAX-RTMPMIN)

! Tables for PARTICLE absorption coefficients

DO IPC=1,N_LAGRANGIAN_CLASSES
   LPC => LAGRANGIAN_PARTICLE_CLASS(IPC)
   IF (LPC%LIQUID_DROPLET) CALL MEAN_CROSS_SECTIONS(PARTICLE_CLASS=IPC)
ENDDO

! Determine angle factors for Lagrangian particles with ORIENTATION
! COSINE_ARRAY holds the cosines of the angles formed by the orientation vector and the radiation directions.
! DLO is the integral of the orientation vector dotted with the directional solid angle of the radiation directions.
! VIEW_ANGLE_FACTOR is the reduction of the radiation due to a view angle less than 180, like a narrow field of view radiometer.

IF (SOLID_PARTICLES) THEN
   ALLOCATE(COSINE_ARRAY(1:NRA))
   ALLOCATE(NEAREST_RADIATION_ANGLE(N_ORIENTATION_VECTOR))
   ALLOCATE(VIEW_ANGLE_FACTOR(N_ORIENTATION_VECTOR))
   VIEW_ANGLE_FACTOR = 0._EB
   DO IO=1,N_ORIENTATION_VECTOR
      DLO = 0._EB
      DO N=1,NRA
         COSINE_ARRAY(N) = ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + &
                           ORIENTATION_VECTOR(2,IO)*DLANG(2,N) + &
                           ORIENTATION_VECTOR(3,IO)*DLANG(3,N)
         IF (-COSINE_ARRAY(N) > COS_HALF_VIEW_ANGLE(IO)) &
            DLO = DLO - (ORIENTATION_VECTOR(1,IO)*DLX(N) + ORIENTATION_VECTOR(2,IO)*DLY(N) + ORIENTATION_VECTOR(3,IO)*DLZ(N))
      ENDDO
      NEAREST_RADIATION_ANGLE(IO) = MINLOC(COSINE_ARRAY,DIM=1)
      VIEW_ANGLE_FACTOR(IO) = PI/DLO
   ENDDO
   DEALLOCATE(COSINE_ARRAY)
ENDIF

! Allocate array needed by angle-specific RADF output files

DO N=1,N_RADF
   RF => RAD_FILE(N)
   ALLOCATE(RF%IL_SAVE(RF%I1:RF%I2,RF%J1:RF%J2,RF%K1:RF%K2,NUMBER_RADIATION_ANGLES))
ENDDO

END SUBROUTINE INIT_RADIATION


!> \brief Compute radiative source term and transfer.
!> \param T Current time (s)
!> \param NM Mesh number
!> \param RAD_ITER Counter for repeated calls to the radiation solver

SUBROUTINE COMPUTE_RADIATION(T,NM,RAD_ITER)

USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP
USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
USE COMPLEX_GEOMETRY
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM,RAD_ITER
REAL(EB), POINTER, DIMENSION(:,:,:) :: KFST4_GAS
REAL(EB) :: TNOW

TNOW=CURRENT_TIME()

CALL POINT_TO_MESH(NM)

IF (RADIATION) THEN
   RADIATION_COMPLETED(NM) = .FALSE.
   CALL RADIATION_FVM
ELSE
   RADIATION_COMPLETED(NM) = .TRUE.
   IF (N_REACTIONS>0) QR = -CHI_R*Q
   IF (INIT_HRRPUV) CALL ADD_VOLUMETRIC_HEAT_SOURCE(0) ! Even if no radiation, add HRRPUV specified on an INIT line
ENDIF

T_USED(9)=T_USED(9)+CURRENT_TIME()-TNOW

CONTAINS


!> \brief Finite Volume Method of radiation transport

SUBROUTINE RADIATION_FVM

USE MIEV
USE MATH_FUNCTIONS, ONLY : INTERPOLATE1D
USE TRAN, ONLY : GET_IJK
USE COMPLEX_GEOMETRY, ONLY : CC_CGSC,CC_SOLID
USE PHYSICAL_FUNCTIONS, ONLY : GET_VOLUME_FRACTION, GET_MASS_FRACTION
REAL(EB) :: RAP, AX, AXU, AXD, AY, AYU, AYD, AZ, AZU, AZD, VC, RU, RD, RP, AFD, &
            ILXU, ILYU, ILZU, QVAL, BBF, BBFA, NCSDROP, RSA_RAT,EFLUX,SOOT_MASS_FRACTION, &
            AIU_SUM,A_SUM,VOL,VC1,AY1,AZ1,DLO,COS_DLO,AILFU, &
            RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC

INTEGER  :: N,NN,IIG,JJG,KKG,I,J,K,IW,ICF,II,JJ,KK,IOR,IC,IWUP,IWDOWN, &
            ISTART, IEND, ISTEP, JSTART, JEND, JSTEP, &
            KSTART, KEND, KSTEP, NSTART, NEND, NSTEP, &
            I_UIID, N_UPDATES, IBND, NOM, ARRAY_INDEX,NRA, &
            IMIN, JMIN, KMIN, IMAX, JMAX, KMAX, N_SLICE, M_IJK, IJK, LL
INTEGER  :: IADD,IFACE,INDCF
INTEGER, ALLOCATABLE :: IJK_SLICE(:,:)
REAL(EB) :: XID,YJD,ZKD,KAPPA_PART_SINGLE,DLF,DLA(3),TSI,TMP_EXTERIOR,TEMP_ORIENTATION(3)
REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_GET
INTEGER :: IID,JJD,KKD,IP
LOGICAL :: UPDATE_INTENSITY
REAL(EB), POINTER, DIMENSION(:,:,:) :: IL,UIIOLD,KAPPA_PART,KFST4_PART,EXTCOE,SCAEFF,SCAEFF_G,IL_UP
REAL(EB), POINTER, DIMENSION(:)     :: OUTRAD_W,INRAD_W,OUTRAD_F,INRAD_F,IL_F
TYPE (OMESH_TYPE), POINTER :: M2
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(VENTS_TYPE), POINTER :: VT
TYPE(RAD_FILE_TYPE), POINTER :: RF
TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE(WALL_TYPE), POINTER :: WC
TYPE(CC_CUTFACE_TYPE), POINTER :: CF
TYPE(CFACE_TYPE), POINTER :: CFA
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_RADIA_TYPE), POINTER :: BR,BR_UP,BR_DOWN
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(INITIALIZATION_TYPE), POINTER :: IN
CHARACTER(20) :: FORMT

! Variables added for the WSGG model
REAL(EB) :: X_H2O, X_CO2, MOL_RAT,PARTIAL_P,R_MIXTURE,TOTAL_P
REAL(EB), ALLOCATABLE, DIMENSION(:) :: Z_ARRAY

ALLOCATE(Z_ARRAY(N_TRACKED_SPECIES))

ALLOCATE( IJK_SLICE(3, IBAR*KBAR) )

KFST4_GAS  => WORK1
IL         => WORK2
UIIOLD     => WORK3
EXTCOE     => WORK4
KAPPA_PART => WORK5
SCAEFF     => WORK6
KFST4_PART => WORK7
IL_UP      => WORK8
SCAEFF_G   => WORK9
OUTRAD_W   => WALL_WORK1
INRAD_W    => WALL_WORK2
IF (CC_IBM) THEN
   OUTRAD_F => FACE_WORK1
   INRAD_F  => FACE_WORK2; INRAD_F  = 0._EB
   IL_F     => FACE_WORK3
ENDIF

! Ratio of solid angle, used in scattering

NRA     = NUMBER_RADIATION_ANGLES
RSA_RAT = 1._EB/(1._EB-1._EB/NRA)

! Check if it time to update radiation intensity field

IF ( MOD(RAD_CALL_COUNTER,TIME_STEP_INCREMENT)==0 .OR. INITIALIZATION_PHASE .OR. ICYC==1 .OR. UPDATE_ALL_ANGLES) THEN
   UPDATE_INTENSITY   = .TRUE.
   EXCHANGE_RADIATION = .TRUE.
ELSE
   UPDATE_INTENSITY   = .FALSE.
   EXCHANGE_RADIATION = .FALSE.
ENDIF

! If this is the last iteration of the solver in a single time step, update the call counter.

IF (RAD_ITER==RADIATION_ITERATIONS) RAD_CALL_COUNTER  = RAD_CALL_COUNTER + 1

! If there is a user-specified HRRPUV on an INIT line, and there are multiple iterations of the radiation solver,
! subtract off the HRRPUV from Q, leaving just the chemical HRR.

IF (INIT_HRRPUV .AND. RAD_ITER>1) CALL ADD_VOLUMETRIC_HEAT_SOURCE(2)

! Initialize the radiative loss to zero for special case models that loop over wavelength bands

IF (WIDE_BAND_MODEL .OR. WSGG_MODEL) THEN
   QR = 0._EB
   IF (STORE_RADIATION_TERMS) THEN
      RADIATION_EMISSION   = 0._EB
      RADIATION_ABSORPTION = 0._EB
   ENDIF
   IF (N_LP_ARRAY_INDICES>0 .AND. PARTICLES_EXISTED) QR_W = 0._EB
ENDIF

! Zero out radiation flux to wall, particles, facets if the intensity is to be updated

IF (UPDATE_INTENSITY) THEN
   DO IW=1,N_INTERNAL_WALL_CELLS+N_EXTERNAL_WALL_CELLS
      WC => WALL(IW)
      IF (WC%B1_INDEX==0 .OR. WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE
      B1 => BOUNDARY_PROP1(WC%B1_INDEX)
      B1%Q_RAD_IN = 0._EB
      SF  => SURFACE(WALL(IW)%SURF_INDEX)
   ENDDO
   DO IP=1,NLP
      LP => LAGRANGIAN_PARTICLE(IP)
      IF (LP%B1_INDEX==0) CYCLE
      B1 => BOUNDARY_PROP1(LP%B1_INDEX)
      B1%Q_RAD_IN = 0._EB
   ENDDO
   DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
      CFA => CFACE(ICF)
      IF (CFA%B1_INDEX==0) CYCLE
      B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
      B1%Q_RAD_IN = 0._EB
      SF  => SURFACE(CFA%SURF_INDEX)
   ENDDO
ENDIF

! Loop over spectral bands

BAND_LOOP: DO IBND = 1,NUMBER_SPECTRAL_BANDS

   KAPPA_PART = 0._EB
   KFST4_GAS  = 0._EB
   KFST4_PART = 0._EB
   SCAEFF = 0._EB
   SCAEFF_G = 0._EB

   ! Calculate fraction on ambient black body radiation

   IF (NUMBER_SPECTRAL_BANDS==1) THEN
      BBFA = 1._EB
   ELSEIF (WSGG_MODEL) THEN
      ! Computing the temperature coefficient in the WSGG model at ambient temperature
      Z_ARRAY(1:N_TRACKED_SPECIES) = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0     ! Mass frac of the tracked species in ambient
      R_MIXTURE = RSUM0                                                           ! Specific gas constant of ambient
      MOL_RAT = GET_VOLUME_FRACTION(H2O_INDEX,Z_ARRAY,R_MIXTURE)/&
         (GET_VOLUME_FRACTION(CO2_INDEX,Z_ARRAY,R_MIXTURE)+TWO_EPSILON_EB) ! Molar ratio
      BBFA = A_WSGG(TMPA,MOL_RAT,IBND)
   ELSE
      BBFA = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMPA)
   ENDIF

   ! Generate water absorption and scattering coefficients

   IF_PARTICLES_INCLUDED: IF (N_LP_ARRAY_INDICES>0) THEN

      IF (NUMBER_SPECTRAL_BANDS==1) THEN
         BBF = 1._EB
         QR_W = 0._EB
      ELSEIF (WSGG_MODEL) THEN
         ! Computing the temperature coefficient in the WSGG model at ambient temperature
         Z_ARRAY(1:N_TRACKED_SPECIES) = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0     ! Mass frac of the tracked species in ambient
         R_MIXTURE = RSUM0                                                           ! Specific gas constant of ambient
         MOL_RAT = GET_VOLUME_FRACTION(H2O_INDEX,Z_ARRAY,R_MIXTURE)/&
            (GET_VOLUME_FRACTION(CO2_INDEX,Z_ARRAY,R_MIXTURE)+TWO_EPSILON_EB) ! Molar ratio
         BBF = A_WSGG(RADTMP,MOL_RAT,IBND)
      ELSE
         BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),RADTMP)
      ENDIF
      IF (LIQUID_DROPLETS .AND. PARTICLES_EXISTED) THEN
         PC_LOOP: DO N=1,N_LAGRANGIAN_CLASSES
            LPC => LAGRANGIAN_PARTICLE_CLASS(N)
            IF (.NOT.LPC%LIQUID_DROPLET) CYCLE PC_LOOP
            ARRAY_INDEX = LPC%ARRAY_INDEX
            IF (ARRAY_INDEX==0) CYCLE PC_LOOP
            DO K=1,KBAR
               DO J=1,JBAR
                  DO I=1,IBAR
                     IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
                     IF (ABS(AVG_DROP_AREA(I,J,K,ARRAY_INDEX))<TWO_EPSILON_EB) CYCLE
                     NCSDROP = AVG_DROP_AREA(I,J,K,ARRAY_INDEX)
                     CALL INTERPOLATE1D(LPC%R50,LPC%WQABS(:,IBND),AVG_DROP_RAD(I,J,K,ARRAY_INDEX),QVAL)
                     KAPPA_PART(I,J,K) = KAPPA_PART(I,J,K) + NCSDROP*QVAL
                     KFST4_PART(I,J,K) = KFST4_PART(I,J,K)+ BBF*NCSDROP*QVAL*FOUR_SIGMA*AVG_DROP_TMP(I,J,K,ARRAY_INDEX)**4
                     CALL INTERPOLATE1D(LPC%R50,LPC%WQSCA(:,IBND),AVG_DROP_RAD(I,J,K,ARRAY_INDEX),QVAL)
                     SCAEFF(I,J,K) = SCAEFF(I,J,K) + NCSDROP*QVAL
                  ENDDO
               ENDDO
            ENDDO
         ENDDO PC_LOOP
      ENDIF

   ENDIF IF_PARTICLES_INCLUDED

   ! Compute the absorption coefficient, KAPPA_PART, for a collection of solid particles

   IF (NLP>0 .AND. SOLID_PARTICLES) THEN
      DO IP = 1,NLP
         LP => LAGRANGIAN_PARTICLE(IP)
         LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
         IF (.NOT.LPC%SOLID_PARTICLE) CYCLE
         B1 => BOUNDARY_PROP1(LP%B1_INDEX)
         BC => BOUNDARY_COORD(LP%BC_INDEX)
         CALL GET_IJK(BC%X,BC%Y,BC%Z,NM,XID,YJD,ZKD,IID,JJD,KKD)
         KAPPA_PART_SINGLE = 0.25_EB*LP%PWT*B1%AREA*LP%RVC*B1%EMISSIVITY
         KAPPA_PART(IID,JJD,KKD) = KAPPA_PART(IID,JJD,KKD) + KAPPA_PART_SINGLE
         KFST4_PART(IID,JJD,KKD) = KFST4_PART(IID,JJD,KKD) + BBF*KAPPA_PART_SINGLE*FOUR_SIGMA*B1%TMP_F**4
      ENDDO
   ENDIF

   ! Compute absorption coefficient of the gases, KAPPA_GAS

   IF (KAPPA0>=0._EB) THEN

      KAPPA_GAS = KAPPA0

   ELSEIF (N_RADCAL_ARRAY_SIZE>0) THEN

      !$OMP PARALLEL PRIVATE(ZZ_GET)
      ALLOCATE(ZZ_GET(1:N_TRACKED_SPECIES))
      !$OMP DO SCHEDULE(STATIC)
      DO K=1,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
               IF (CC_IBM) THEN
                  IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
               ENDIF
               ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)
               KAPPA_GAS(I,J,K) = GET_KAPPA(ZZ_GET,TMP(I,J,K),IBND)
            ENDDO
         ENDDO
      ENDDO
      !$OMP END DO
      DEALLOCATE(ZZ_GET)
      !$OMP END PARALLEL
   ENDIF

   ! Compute source term KAPPA*4*SIGMA*TMP**4

   WIDE_BAND_MODEL_IF: IF (WIDE_BAND_MODEL) THEN

      ! Wide band model

      DO K=1,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
               IF (CC_IBM) THEN
                  IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
               ENDIF
               BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMP(I,J,K))
               KFST4_GAS(I,J,K) = BBF*KAPPA_GAS(I,J,K)*FOUR_SIGMA*TMP(I,J,K)**4
            ENDDO
         ENDDO
      ENDDO

   ELSEIF (WSGG_MODEL) THEN WIDE_BAND_MODEL_IF

      ALPHA_CC = 1._EB
      DO K=1,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
               IF (CC_IBM) THEN
                  ALPHA_CC = 1._EB
                  IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
                  IC = CCVAR(I,J,K,CC_IDCC)
                  IF (IC>0) ALPHA_CC = CUT_CELL(IC)%ALPHA_CC
               ENDIF
               Z_ARRAY(1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)                  ! Mass frac of the tracked species
               R_MIXTURE = RSUM(I,J,K)                                                       ! Specific gas constant of the mixture
               X_H2O = GET_VOLUME_FRACTION(H2O_INDEX,Z_ARRAY,R_MIXTURE)
               X_CO2 = GET_VOLUME_FRACTION(CO2_INDEX,Z_ARRAY,R_MIXTURE)
               MOL_RAT = X_H2O/(X_CO2 + TWO_EPSILON_EB)                                      ! Molar ratio
               TOTAL_P = PBAR(K,PRESSURE_ZONE(I,J,K)) + RHO(I,J,K)*(H(I,J,K)-KRES(I,J,K))    ! Total pressure
               PARTIAL_P = TOTAL_P*(X_H2O + X_CO2)/P_STP                                     ! Partial press of the CO2-H2O mixture
               BBF = A_WSGG(TMP(I,J,K),MOL_RAT,IBND)                                         ! Temp coefficient for the jth gas
               CALL GET_MASS_FRACTION(Z_ARRAY,SOOT_INDEX,SOOT_MASS_FRACTION)

               KAPPA_GAS(I,J,K) = KAPPA_WSGG(X_H2O, X_CO2,MOL_RAT,PARTIAL_P,IBND) + &
                                  KAPPA_SOOT(RHO(I,J,K)*SOOT_MASS_FRACTION,TMP(I,J,K))       ! Absorp coeff for the jth gas

               KFST4_GAS(I,J,K) = BBF*KAPPA_GAS(I,J,K)*FOUR_SIGMA*TMP(I,J,K)**4._EB
               IF (CHI_R(I,J,K)*Q(I,J,K)>QR_CLIP) THEN ! Precomputation of quantities for the RTE source term correction
                     VOL = R(I)*DX(I)*DY(J)*DZ(K)*ALPHA_CC
                     RAD_Q_SUM = RAD_Q_SUM + (BBF*CHI_R(I,J,K)*Q(I,J,K) + KAPPA_GAS(I,J,K)*UIID(I,J,K,IBND))*VOL
                     KFST4_SUM = KFST4_SUM + KFST4_GAS(I,J,K)*VOL
               ENDIF
            ENDDO
         ENDDO
      ENDDO

      ! Correct the source term in the RTE based on user-specified RADIATIVE_FRACTION on REAC

      DO K=1,KBAR
         DO J=1,JBAR
            DO I=1,IBAR
               IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
               IF (CC_IBM) THEN
                  IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
               ENDIF
               IF (CHI_R(I,J,K)*Q(I,J,K)>QR_CLIP) KFST4_GAS(I,J,K) = KFST4_GAS(I,J,K)*RTE_SOURCE_CORRECTION_FACTOR
            ENDDO
         ENDDO
      ENDDO

   ELSE WIDE_BAND_MODEL_IF

      ! Gray gas model

      RTE_SOURCE_CORRECTION_IF: IF (RTE_SOURCE_CORRECTION) THEN ! default RTE_SOURCE_CORRECTION=.TRUE.

         ! Only apply the correction to KFST4_GAS for gray gas model

         !$OMP PARALLEL PRIVATE(RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC) SHARED(RAD_Q_SUM,KFST4_SUM)

         RAD_Q_SUM_PARTIAL = 0._EB
         KFST4_SUM_PARTIAL = 0._EB
         ALPHA_CC = 1._EB

         !$OMP DO PRIVATE(VOL,IC)
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
                  IF (CC_IBM) THEN
                     ALPHA_CC = 1._EB
                     IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
                     IC = CCVAR(I,J,K,CC_IDCC)
                     IF (IC>0) ALPHA_CC = CUT_CELL(IC)%ALPHA_CC
                  ENDIF
                  KFST4_GAS(I,J,K) = KAPPA_GAS(I,J,K)*FOUR_SIGMA*TMP(I,J,K)**4
                  IF (CHI_R(I,J,K)*Q(I,J,K)>QR_CLIP) THEN
                     VOL = R(I)*DX(I)*DY(J)*DZ(K)*ALPHA_CC
                     RAD_Q_SUM_PARTIAL = RAD_Q_SUM_PARTIAL + (CHI_R(I,J,K)*Q(I,J,K)+KAPPA_GAS(I,J,K)*UII(I,J,K))*VOL
                     KFST4_SUM_PARTIAL = KFST4_SUM_PARTIAL + KFST4_GAS(I,J,K)*VOL
                  ENDIF
               ENDDO
            ENDDO
         ENDDO
         !$OMP END DO

         !$OMP CRITICAL
         RAD_Q_SUM = RAD_Q_SUM + RAD_Q_SUM_PARTIAL
         KFST4_SUM = KFST4_SUM + KFST4_SUM_PARTIAL
         !$OMP END CRITICAL

         ! Correct the source term in the RTE based on user-specified RADIATIVE_FRACTION on REAC

         !$OMP DO
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
                  IF (CC_IBM) THEN
                     IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
                  ENDIF
                  IF (CHI_R(I,J,K)*Q(I,J,K)>QR_CLIP) KFST4_GAS(I,J,K) = KFST4_GAS(I,J,K)*RTE_SOURCE_CORRECTION_FACTOR
               ENDDO
            ENDDO
         ENDDO
         !$OMP END DO

         !$OMP END PARALLEL

      ELSE RTE_SOURCE_CORRECTION_IF  ! OPTICALLY_THIN

         ! Use specified radiative fraction

         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
                  IF (CC_IBM) THEN
                     IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
                  ENDIF
                  KFST4_GAS(I,J,K) = CHI_R(I,J,K)*Q(I,J,K)+KAPPA_GAS(I,J,K)*UII(I,J,K)
               ENDDO
            ENDDO
         ENDDO

      ENDIF RTE_SOURCE_CORRECTION_IF

   ENDIF WIDE_BAND_MODEL_IF

   ! Add contribution to source term from a user-specified volumetric heat release rate

   IF (INIT_HRRPUV) CALL ADD_VOLUMETRIC_HEAT_SOURCE(1)

   ! Compute the added contribution of any condensed species

   IF (ANY(SPECIES_MIXTURE%EVAPORATION_SMIX_INDEX > 0)) THEN

      IF (NUMBER_SPECTRAL_BANDS==1) THEN
         BBF = 1._EB
      ELSE
         BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),RADTMP)
      ENDIF

      SLOOP: DO N = 1, N_TRACKED_SPECIES
         IF (SPECIES_MIXTURE(N)%EVAPORATION_SMIX_INDEX <= 0) CYCLE SLOOP
         DO K=1,KBAR
            DO J=1,JBAR
               DO I=1,IBAR
                  IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
                  IF (CC_IBM) THEN
                     IF (CCVAR(I,J,K,CC_CGSC)==CC_SOLID) CYCLE
                  ENDIF
                  IF (ZZ(I,J,K,N) < TWO_EPSILON_EB) CYCLE
                  NCSDROP = 1.5_EB*ZZ(I,J,K,N)*RHO(I,J,K)/ &
                            (SPECIES(SPECIES_MIXTURE(N)%SINGLE_SPEC_INDEX)%DENSITY_LIQUID*SPECIES_MIXTURE(N)%MEAN_DIAMETER)
                  CALL INTERPOLATE1D(SPECIES_MIXTURE(N)%R50,SPECIES_MIXTURE(N)%WQABS(:,IBND),&
                                     0.5_EB*SPECIES_MIXTURE(N)%MEAN_DIAMETER,QVAL)
                  KAPPA_GAS(I,J,K) = KAPPA_GAS(I,J,K) + NCSDROP*QVAL
                  KFST4_GAS(I,J,K) = KFST4_GAS(I,J,K) + BBF*NCSDROP*QVAL*FOUR_SIGMA*TMP(I,J,K)**4
                  CALL INTERPOLATE1D(SPECIES_MIXTURE(N)%R50,SPECIES_MIXTURE(N)%WQSCA(:,IBND),&
                                     0.5_EB*SPECIES_MIXTURE(N)%MEAN_DIAMETER,QVAL)
                  SCAEFF_G(I,J,K) = SCAEFF_G(I,J,K) + NCSDROP*QVAL
               ENDDO
            ENDDO
         ENDDO
      ENDDO SLOOP
   ENDIF

   ! Calculate extinction coefficient

   EXTCOE = KAPPA_GAS + KAPPA_PART + (SCAEFF+SCAEFF_G)*RSA_RAT

   ! Update intensity field

   INTENSITY_UPDATE: IF (UPDATE_INTENSITY) THEN

      IF (WIDE_BAND_MODEL .OR. WSGG_MODEL) THEN
         UIIOLD = UIID(:,:,:,IBND)
      ELSE
         UIIOLD = UII
      ENDIF
      UII = 0._EB

      ! Compute boundary condition intensity emissivity*sigma*Tw**4/pi or emissivity*QRADOUT/pi for wall with internal radiation

      BBF = 1.0_EB
      DO IW = 1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
         WC => WALL(IW)
         IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE
         B1 => BOUNDARY_PROP1(WC%B1_INDEX)
         BC    => BOUNDARY_COORD(WC%BC_INDEX)
         IF (WC%BOUNDARY_TYPE == OPEN_BOUNDARY) THEN
            OUTRAD_W(IW) = BBFA*RPI*SIGMA*TMPA**4
         ELSE
            IF (WIDE_BAND_MODEL) BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),B1%TMP_F)
            IF (WSGG_MODEL) THEN
               Z_ARRAY(1:N_TRACKED_SPECIES) = ZZ(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES)
               R_MIXTURE = RSUM(BC%IIG,BC%JJG,BC%KKG)
               MOL_RAT = GET_VOLUME_FRACTION(H2O_INDEX,Z_ARRAY,R_MIXTURE)/&
                  (GET_VOLUME_FRACTION(CO2_INDEX,Z_ARRAY,R_MIXTURE) + TWO_EPSILON_EB)
               BBF = A_WSGG(B1%TMP_F,MOL_RAT,IBND) ! Temperature coefficient for the jth gray gas in the boundary
            ENDIF                                     ! (use information of the cell adjacent to the boundary)
            SF  => SURFACE(WC%SURF_INDEX)
            OUTRAD_W(IW) = BBF*RPI*B1%Q_RAD_OUT
         ENDIF
      ENDDO

      BBF = 1.0_EB
      DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
         CFA => CFACE(ICF)
         B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
         BC    => BOUNDARY_COORD(CFA%BC_INDEX)
         IF (WIDE_BAND_MODEL) BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),B1%TMP_F)
         SF => SURFACE(CFA%SURF_INDEX)
         OUTRAD_F(ICF) = BBF*RPI*B1%Q_RAD_OUT
      ENDDO

      ! Compute boundary condition term incoming radiation integral

      DO IW = 1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
         WC => WALL(IW)
         IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE
         BC => BOUNDARY_COORD(WC%BC_INDEX)
         BR => BOUNDARY_RADIA(WC%BR_INDEX)
         INRAD_W(IW) = SUM(-DLN(BC%IOR,:)*BR%BAND(IBND)%ILW(:), 1, DLN(BC%IOR,:)<0._EB)
      ENDDO

      DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
         CFA => CFACE(ICF)
         BR  => BOUNDARY_RADIA(CFA%BR_INDEX)
         BC  => BOUNDARY_COORD(CFA%BC_INDEX)
         DO N=1,NRA
            DLA = (/DLX(N),DLY(N),DLZ(N)/)
            DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle
            IF (DLF<0._EB) INRAD_F(ICF) = INRAD_F(ICF) - DLF*BR%BAND(IBND)%ILW(N)
         ENDDO
      ENDDO

      ! If updating intensities first time, sweep ALL angles

      N_UPDATES = 1
      IF (INITIALIZATION_PHASE .OR. ICYC==1 .OR. UPDATE_ALL_ANGLES) N_UPDATES = ANGLE_INCREMENT

      UPDATE_LOOP: DO I_UIID = 1,N_UPDATES

         ! Update counters inside the radiation routine

         ANGLE_INC_COUNTER = MOD(ANGLE_INC_COUNTER,ANGLE_INCREMENT) + 1

         ! If this is the last set of angles to update, indicate that the radiation routine has finished a full update

         IF (ANGLE_INC_COUNTER==ANGLE_INCREMENT) RADIATION_COMPLETED(NM) = .TRUE.

         ! Zero out UIID, the integrated intensity

         IF (WIDE_BAND_MODEL .OR. WSGG_MODEL) THEN
            UIID(:,:,:,IBND) = 0._EB
         ELSE
            UIID(:,:,:,ANGLE_INC_COUNTER) = 0._EB
         ENDIF

         DO IW=1,N_EXTERNAL_WALL_CELLS
            WC => WALL(IW)
            IF (WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) CYCLE
            BR => BOUNDARY_RADIA(WC%BR_INDEX)
            BR%BAND(IBND)%ILW(ANGLE_INC_COUNTER) = 0._EB
         ENDDO

         ! Set the bounds and increment for the angleloop. Step downdard because in cylindrical case the Nth angle
         ! boundary condition comes from (N+1)th angle.

         NSTART    = NRA - ANGLE_INC_COUNTER + 1
         NEND      = 1
         NSTEP     = -ANGLE_INCREMENT

         IL(:,:,:) = BBFA*RPI_SIGMA*TMPA4

         ANGLE_LOOP: DO N = NSTART,NEND,NSTEP  ! Sweep through control angles

            CELL_ILW = -HUGE(EB)

            ! Boundary conditions: Intensities leaving the boundaries.

            !$OMP PARALLEL DO PRIVATE(WC,BC,BR,B1,IOR,II,JJ,KK,LL,NOM,VT,TSI,TMP_EXTERIOR,IC) SCHEDULE(GUIDED)
            WALL_LOOP1: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
               WC => WALL(IW)
               IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP1
               BC => BOUNDARY_COORD(WC%BC_INDEX)
               BR => BOUNDARY_RADIA(WC%BR_INDEX)
               B1 => BOUNDARY_PROP1(WC%B1_INDEX)
               IOR = BC%IOR
               IF (DLN(IOR,N) < 0._EB) CYCLE WALL_LOOP1
               II  = BC%II
               JJ  = BC%JJ
               KK  = BC%KK
               IF (.NOT.TWO_D .OR. ABS(IOR)/=2) THEN
                  SELECT CASE (WC%BOUNDARY_TYPE)
                     CASE (OPEN_BOUNDARY)
                        VT => VENTS(WC%VENT_INDEX)
                        IF (VT%TMP_EXTERIOR>0._EB) THEN
                           TSI = T - T_BEGIN
                           TMP_EXTERIOR = TMP_0(KK)+EVALUATE_RAMP(TSI,VT%TMP_EXTERIOR_RAMP_INDEX)*(VT%TMP_EXTERIOR-TMP_0(KK))
                           IL(II,JJ,KK) = BBFA*RPI_SIGMA*TMP_EXTERIOR**4
                        ELSE
                           IL(II,JJ,KK) = BBFA*RPI_SIGMA*TMPA4
                        ENDIF
                     CASE (MIRROR_BOUNDARY)
                        BR%BAND(IBND)%ILW(N) = BR%BAND(IBND)%ILW(DLM(N,ABS(IOR)))
                        IL(II,JJ,KK) = BR%BAND(IBND)%ILW(N)
                     CASE (INTERPOLATED_BOUNDARY)
                        ! IL_R holds the intensities from mesh NOM in the ghost cells of mesh NM.
                        ! IL(II,JJ,KK) is the average of the intensities from the other mesh.
                        NOM = EXTERNAL_WALL(IW)%NOM
                        IL(II,JJ,KK) = 0._EB
                        DO LL=EXTERNAL_WALL(IW)%NIC_MIN,EXTERNAL_WALL(IW)%NIC_MAX
                           IL(II,JJ,KK) = IL(II,JJ,KK) + OMESH(NOM)%IL_R(LL,N,IBND)
                        ENDDO
                        IL(II,JJ,KK) = IL(II,JJ,KK)/REAL(EXTERNAL_WALL(IW)%NIC_MAX-EXTERNAL_WALL(IW)%NIC_MIN+1,EB)
                     CASE DEFAULT ! solid wall
                        BR%BAND(IBND)%ILW(N) = OUTRAD_W(IW) + RPI*(1._EB-B1%EMISSIVITY)*INRAD_W(IW)
                        IC = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG)
                        CELL_ILW(IC,ABS(IOR)) = BR%BAND(IBND)%ILW(N)
                  END SELECT
               ELSEIF (CYLINDRICAL) THEN
                  IF (WC%BOUNDARY_TYPE==OPEN_BOUNDARY) CYCLE WALL_LOOP1
                  IL(II,JJ,KK) = BR%BAND(IBND)%ILW(N)
               ENDIF
            ENDDO WALL_LOOP1
            !$OMP END PARALLEL DO

            DLA = (/DLX(N),DLY(N),DLZ(N)/)
            CFACE_LOOP1: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
               CFA => CFACE(ICF)
               IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP1
               BR  => BOUNDARY_RADIA(CFA%BR_INDEX)
               BC => BOUNDARY_COORD(CFA%BC_INDEX)
               B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
               DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle
               IF (DLF<0._EB) CYCLE CFACE_LOOP1
               BR%BAND(IBND)%ILW(N) = OUTRAD_F(ICF) + RPI*(1._EB-B1%EMISSIVITY)*INRAD_F(ICF)
            ENDDO CFACE_LOOP1

            ! Determine sweep direction in physical space

            ISTART = 1
            JSTART = 1
            KSTART = 1
            IEND   = IBAR
            JEND   = JBAR
            KEND   = KBAR
            ISTEP  = 1
            JSTEP  = 1
            KSTEP  = 1
            IMIN = ISTART
            JMIN = JSTART
            KMIN = KSTART
            IMAX = IEND
            JMAX = JEND
            KMAX = KEND
            IF (DLX(N) < 0._EB) THEN
               ISTART = IBAR
               IEND   = 1
               ISTEP  = -1
               IMIN = IEND
               IMAX = ISTART
            ENDIF
            IF (DLY(N) < 0._EB) THEN
               JSTART = JBAR
               JEND   = 1
               JSTEP  = -1
               JMIN = JEND
               JMAX = JSTART
            ENDIF
            IF (DLZ(N) < 0._EB) THEN
               KSTART = KBAR
               KEND   = 1
               KSTEP  = -1
               KMIN = KEND
               KMAX = KSTART
            ENDIF

            GEOMETRY2: IF (CYLINDRICAL) THEN  ! Sweep in axisymmetric geometry
               J = 1
               CKLOOP: DO K=KSTART,KEND,KSTEP
                  CILOOP: DO I=ISTART,IEND,ISTEP
                     ILXU = IL(I-ISTEP,J,K)
                     ILYU = IL(I,J-JSTEP,K)
                     ILZU = IL(I,J,K-KSTEP)
                     IF (DLX(N)>=0._EB) THEN
                        RU  = R(I-1)
                        RD  = R(I)
                     ELSE
                        RU  = R(I)
                        RD  = R(I-1)
                     ENDIF
                     RP  = SQRT(0.5_EB*(RU**2+RD**2))
                     VC  = DX(I)  * RP*DPHI0 * DZ(K)
                     AXU = 2._EB*SIN(DPHI0/2.)*RU       * DZ(K) * ABS(DLX(N))
                     AXD = 2._EB*SIN(DPHI0/2.)*RD       * DZ(K) * ABS(DLX(N))
                     AYU = DX(I)             * DZ(K) * ABS(DLB(N))
                     AYD = DX(I)             * DZ(K) * ABS(DLY(N))
                     AZ  = DX(I)  * RP*DPHI0         * ABS(DLZ(N))
                     IF (MODULO(N,NRP(1))==1) AYD = 0._EB  ! Zero out the terms involving symmetric overhang
                     IF (MODULO(N,NRP(1))==0) AYU = 0._EB
                     IC = CELL_INDEX(I,J,K)
                     IF (IC/=0) THEN
                        IF (CELL(IC)%SOLID) CYCLE CILOOP
                        IF (CELL_ILW(IC,1)>-1.E6_EB) ILXU = CELL_ILW(IC,1)
                        IF (CELL_ILW(IC,2)>-1.E6_EB) ILYU = CELL_ILW(IC,2)
                        IF (CELL_ILW(IC,3)>-1.E6_EB) ILZU = CELL_ILW(IC,3)
                     ENDIF
                     AIU_SUM = AXU*ILXU + AYU*ILYU + AZ*ILZU
                     A_SUM = AXD + AYD + AZ
                     RAP = 1._EB/(A_SUM + EXTCOE(I,J,K)*VC*RSA(N))
                     IL(I,J,K) = MAX(0._EB, RAP * (AIU_SUM + VC*RSA(N)*RFPI* &
                                 ( KFST4_GAS(I,J,K) + KFST4_PART(I,J,K) + RSA_RAT*&
                                 (SCAEFF(I,J,K)+SCAEFF_G(I,J,K))*UIIOLD(I,J,K) ) ) )
                     IF (SOLID_PARTICLES) IL_UP(I,J,K) = MAX(0._EB,AIU_SUM/A_SUM)
                  ENDDO CILOOP
               ENDDO CKLOOP

            ELSEIF (TWO_D) THEN GEOMETRY2  ! Sweep in 2D cartesian geometry
               J = 1
               K2LOOP: DO K=KSTART,KEND,KSTEP
                  I2LOOP: DO I=ISTART,IEND,ISTEP
                     ILXU  = IL(I-ISTEP,J,K)
                     ILZU  = IL(I,J,K-KSTEP)
                     VC  = DX(I) * DZ(K)
                     AX  =         DZ(K) * ABS(DLX(N))
                     AZ  = DX(I)         * ABS(DLZ(N))
                     IC = CELL_INDEX(I,J,K)
                     IF (IC/=0) THEN
                        IF (CELL(IC)%SOLID) CYCLE I2LOOP
                        IF (CELL_ILW(IC,1)>-1.E6_EB) ILXU = CELL_ILW(IC,1)
                        IF (CELL_ILW(IC,3)>-1.E6_EB) ILZU = CELL_ILW(IC,3)
                     ENDIF
                     AIU_SUM = AX*ILXU + AZ*ILZU
                     A_SUM = AX + AZ
                     RAP = 1._EB/(A_SUM + EXTCOE(I,J,K)*VC*RSA(N))
                     IL(I,J,K) = MAX(0._EB, RAP * (AIU_SUM + VC*RSA(N)*RFPI* &
                                    (KFST4_GAS(I,J,K) + KFST4_PART(I,J,K) + RSA_RAT* &
                                    (SCAEFF(I,J,K)+SCAEFF_G(I,J,K))*UIIOLD(I,J,K) ) ) )
                     IF (SOLID_PARTICLES) IL_UP(I,J,K) = MAX(0._EB,AIU_SUM/A_SUM)
                  ENDDO I2LOOP
               ENDDO K2LOOP

            ELSE GEOMETRY2  ! Sweep in 3D cartesian geometry

               IPROP_LOOP: DO N_SLICE = ISTEP*ISTART + JSTEP*JSTART + KSTEP*KSTART, &
                                        ISTEP*IEND + JSTEP*JEND + KSTEP*KEND
                  M_IJK = 0
                  DO K = KMIN, KMAX
                     IF (ISTEP*JSTEP > 0) THEN ! I STARTS HIGH
                        JSTART = MAX(JMIN, JSTEP*(N_SLICE - KSTEP*K - ISTEP*IMAX))
                        JEND   = MIN(JMAX, JSTEP*(N_SLICE - KSTEP*K - ISTEP*IMIN))
                     ELSE IF (ISTEP*JSTEP < 0) THEN ! I STARTS LOW
                        JSTART = MAX(JMIN, JSTEP*(N_SLICE - KSTEP*K - ISTEP*IMIN))
                        JEND   = MIN(JMAX, JSTEP*(N_SLICE - KSTEP*K - ISTEP*IMAX))
                     ENDIF
                     IF (JSTART > JEND) THEN
                        CYCLE
                     ENDIF
                     DO J = JSTART, JEND
                        I = ISTEP * (N_SLICE - J*JSTEP - K*KSTEP)
                        M_IJK = M_IJK+1
                        IJK_SLICE(:,M_IJK) = (/I,J,K/)
                     ENDDO
                  ENDDO

                  !$OMP PARALLEL DO SCHEDULE(GUIDED) &
                  !$OMP& PRIVATE(I, J, K, AY1, AX, VC1, AZ1, IC, ILXU, ILYU, AILFU, &
                  !$OMP& ILZU, VC, AY, AZ, AXU, AYU, AZU, AXD, AYD, AZD, AFD, &
                  !$OMP& IW, WC, BR, CF, CFA, BC, DLF, A_SUM, AIU_SUM, RAP, &
                  !$OMP& ICF, INDCF, IADD, IFACE )

                  SLICE_LOOP: DO IJK = 1, M_IJK
                     I = IJK_SLICE(1,IJK)
                     J = IJK_SLICE(2,IJK)
                     K = IJK_SLICE(3,IJK)

                     AY1 = DZ(K) * ABS(DLY(N))
                     AX  = DY(J) * DZ(K) * ABS(DLX(N))
                     VC1 = DY(J) * DZ(K)
                     AZ1 = DY(J) * ABS(DLZ(N))
                     ILXU  = IL(I-ISTEP,J,K)
                     ILYU  = IL(I,J-JSTEP,K)
                     ILZU  = IL(I,J,K-KSTEP)
                     IC = CELL_INDEX(I,J,K)
                     IF (IC/=0) THEN
                        IF (CELL(IC)%SOLID) CYCLE SLICE_LOOP
                        IF (CELL_ILW(IC,1)>-1.E6_EB) ILXU = CELL_ILW(IC,1)
                        IF (CELL_ILW(IC,2)>-1.E6_EB) ILYU = CELL_ILW(IC,2)
                        IF (CELL_ILW(IC,3)>-1.E6_EB) ILZU = CELL_ILW(IC,3)
                     ENDIF
                     VC  = DX(I) * VC1
                     AY  = DX(I) * AY1
                     AZ  = DX(I) * AZ1

                     ! With geom, downwind faces need to be separate from upwind
                     AXU = AX; AXD = AX
                     AYU = AY; AYD = AY
                     AZU = AZ; AZD = AZ

                     AILFU = 0._EB; AFD = 0._EB

                     ! Cut face contributions
                     IF (CC_IBM) THEN
                        IF (CCVAR(I,J,K,CC_CGSC) == CC_SOLID) CYCLE SLICE_LOOP
                        INDCF = CCVAR(I,J,K,CC_IDCF)
                        IF (INDCF>0) THEN ! otherwise can assume CC_GAS?
                           CF => CUT_FACE(INDCF)
                           ! Loop through CFACES and assign as upwind or downwind
                           DO IFACE=1,CF%NFACE
                              CFA => CFACE(CF%CFACE_INDEX(IFACE))
                              BC => BOUNDARY_COORD(CFA%BC_INDEX)
                              DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle
                              IF (DLF>0._EB) THEN ! upwind
                                 BR  => BOUNDARY_RADIA(CFA%BR_INDEX)
                                 AILFU = AILFU + DLF*CFA%AREA*BR%BAND(IBND)%ILW(N)
                              ELSEIF (DLF<0._EB) THEN ! downwind
                                 AFD = AFD + ABS(DLF)*CFA%AREA
                              ENDIF
                           ENDDO

                           ! Adjust gas phase areas based on apertures
                           ! X axis
                           IADD= -(1+ISTEP)/2 ! upwind
                           ICF = FCVAR(I+IADD,J,K,CC_IDCF,IAXIS)
                           IF (ICF>0) AXU = AXU*CUT_FACE(ICF)%ALPHA_CF
                           IF(FCVAR(I+IADD,J,K,CC_FGSC,IAXIS)==CC_SOLID) AXU = 0._EB
                           IADD= (ISTEP-1)/2 ! downwind
                           ICF = FCVAR(I+IADD,J,K,CC_IDCF,IAXIS)
                           IF (ICF>0) AXD = AXD*CUT_FACE(ICF)%ALPHA_CF
                           IF(FCVAR(I+IADD,J,K,CC_FGSC,IAXIS)==CC_SOLID) AXD = 0._EB
                           ! Y axis
                           IADD= -(1+JSTEP)/2 ! upwind
                           ICF = FCVAR(I,J+IADD,K,CC_IDCF,JAXIS)
                           IF (ICF>0) AYU = AYU*CUT_FACE(ICF)%ALPHA_CF
                           IF(FCVAR(I,J+IADD,K,CC_FGSC,JAXIS)==CC_SOLID) AYU = 0._EB
                           IADD= (JSTEP-1)/2 ! downwind
                           ICF = FCVAR(I,J+IADD,K,CC_IDCF,JAXIS)
                           IF (ICF>0) AYD = AYD*CUT_FACE(ICF)%ALPHA_CF
                           IF(FCVAR(I,J+IADD,K,CC_FGSC,JAXIS)==CC_SOLID) AYD = 0._EB
                           ! Z axis
                           IADD= -(1+KSTEP)/2 ! upwind
                           ICF = FCVAR(I,J,K+IADD,CC_IDCF,KAXIS)
                           IF (ICF>0) AZU = AZU*CUT_FACE(ICF)%ALPHA_CF
                           IF(FCVAR(I,J,K+IADD,CC_FGSC,KAXIS)==CC_SOLID) AZU = 0._EB
                           IADD= (KSTEP-1)/2 ! downwind
                           ICF = FCVAR(I,J,K+IADD,CC_IDCF,KAXIS)
                           IF (ICF>0) AZD = AZD*CUT_FACE(ICF)%ALPHA_CF
                           IF(FCVAR(I,J,K+IADD,CC_FGSC,KAXIS)==CC_SOLID) AZD = 0._EB

                           ! Adjust volume
                           IC = CCVAR(I,J,K,CC_IDCC)
                           VC = VC*CUT_CELL(IC)%ALPHA_CC

                        ENDIF
                     ENDIF

                     A_SUM = AXD + AYD + AZD + AFD
                     AIU_SUM = AXU*ILXU + AYU*ILYU + AZU*ILZU + AILFU
                     IF (SOLID_PARTICLES) IL_UP(I,J,K) = MAX(0._EB,AIU_SUM/A_SUM)
                     RAP = 1._EB/(A_SUM + EXTCOE(I,J,K)*VC*RSA(N))
                     IL(I,J,K) = MAX(0._EB, RAP * (AIU_SUM + VC*RSA(N)*RFPI* &
                                     ( KFST4_GAS(I,J,K) + KFST4_PART(I,J,K) + RSA_RAT*&
                                     (SCAEFF(I,J,K)+SCAEFF_G(I,J,K))*UIIOLD(I,J,K) ) ) )

                  ENDDO SLICE_LOOP
                  !$OMP END PARALLEL DO

               ENDDO IPROP_LOOP

            ENDIF GEOMETRY2

            ! Copy the Y-downwind intensities to Y-upwind in cylindrical case

            IF (CYLINDRICAL) THEN
               J=1
               CKLOOP2: DO K=1,KBAR
               CILOOP2: DO I=1,IBAR
                  IC = CELL_INDEX(I,J,K)
                  IF (CELL(IC)%SOLID) CYCLE CILOOP2
                  IWUP   = CELL(CELL_INDEX(I,J,K))%WALL_INDEX(-2)
                  IWDOWN = CELL(CELL_INDEX(I,J,K))%WALL_INDEX( 2)
                  IF (IWUP /=0 .AND. IWDOWN /= 0) THEN
                     BR_UP   => BOUNDARY_RADIA(WALL(IWUP)%BR_INDEX)
                     BR_DOWN => BOUNDARY_RADIA(WALL(IWDOWN)%BR_INDEX)
                     IF (MODULO(N,NRP(1))==1) THEN
                        BR_UP%BAND(IBND)%ILW(N)   = BR_DOWN%BAND(IBND)%ILW(N)
                     ELSE
                        BR_UP%BAND(IBND)%ILW(N-1) = BR_DOWN%BAND(IBND)%ILW(N)
                     ENDIF
                  ENDIF
               ENDDO CILOOP2
               ENDDO CKLOOP2
            ENDIF

            ! Boundary values: Incoming radiation

            !$OMP PARALLEL PRIVATE(IOR, IIG, JJG, KKG, WC, BC, BR)
            !$OMP DO SCHEDULE(GUIDED)
            WALL_LOOP2: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
               WC => WALL(IW)
               IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP2
               IF (WC%BOUNDARY_TYPE==OPEN_BOUNDARY) CYCLE WALL_LOOP2
               BC => BOUNDARY_COORD(WC%BC_INDEX)
               BR => BOUNDARY_RADIA(WC%BR_INDEX)
               IOR = BC%IOR
               IF (TWO_D .AND. .NOT.CYLINDRICAL  .AND. ABS(IOR)==2) CYCLE WALL_LOOP2  ! 2-D non cylindrical
               IF (DLN(IOR,N)>=0._EB) CYCLE WALL_LOOP2     ! outgoing
               IIG = BC%IIG
               JJG = BC%JJG
               KKG = BC%KKG
               INRAD_W(IW) = INRAD_W(IW) + DLN(IOR,N) * BR%BAND(IBND)%ILW(N) ! update incoming rad, step 1
               BR%BAND(IBND)%ILW(N) = IL(IIG,JJG,KKG)
               INRAD_W(IW) = INRAD_W(IW) - DLN(IOR,N) * BR%BAND(IBND)%ILW(N) ! update incoming rad, step 2
            ENDDO WALL_LOOP2
            !$OMP END DO

            !$OMP DO SCHEDULE(GUIDED)
            WALL_LOOP3: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
               WC => WALL(IW)
               IF (WC%BOUNDARY_TYPE/=OPEN_BOUNDARY)   CYCLE WALL_LOOP3
               BC => BOUNDARY_COORD(WC%BC_INDEX)
               IF (DLN(BC%IOR,N)>=0._EB) CYCLE WALL_LOOP3     ! outgoing
               BR => BOUNDARY_RADIA(WC%BR_INDEX)
               BR%BAND(IBND)%ILW(ANGLE_INC_COUNTER) = BR%BAND(IBND)%ILW(ANGLE_INC_COUNTER) - DLN(BC%IOR,N)*IL(BC%IIG,BC%JJG,BC%KKG)
            ENDDO WALL_LOOP3
            !$OMP END DO
            !$OMP END PARALLEL

            CFACE_LOOP2: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
               CFA => CFACE(ICF)
               IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP2
               BC => BOUNDARY_COORD(CFA%BC_INDEX)
               DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle
               IF (DLF>=0._EB) CYCLE CFACE_LOOP2      ! outgoing
               BR  => BOUNDARY_RADIA(CFA%BR_INDEX)
               BC => BOUNDARY_COORD(CFA%BC_INDEX)
               INRAD_F(ICF) = INRAD_F(ICF) + DLF * BR%BAND(IBND)%ILW(N) ! update incoming rad, step 1
               BR%BAND(IBND)%ILW(N) = IL(BC%IIG,BC%JJG,BC%KKG)
               INRAD_F(ICF) = INRAD_F(ICF) - DLF * BR%BAND(IBND)%ILW(N) ! update incoming rad, step 2
            ENDDO CFACE_LOOP2
            CFACE_LOOP3: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
               CFA => CFACE(ICF)
               IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP3
               BC => BOUNDARY_COORD(CFA%BC_INDEX)
               DLF = DOT_PRODUCT(BC%NVEC,DLA) ! face normal * radiation angle
               IF (DLF>=0._EB) CYCLE CFACE_LOOP3      ! outgoing
               BR  => BOUNDARY_RADIA(CFA%BR_INDEX)
               BC => BOUNDARY_COORD(CFA%BC_INDEX)
               BR%BAND(IBND)%ILW(ANGLE_INC_COUNTER) = BR%BAND(IBND)%ILW(ANGLE_INC_COUNTER) - DLF*IL(BC%IIG,BC%JJG,BC%KKG)
            ENDDO CFACE_LOOP3
            ! ENDIF

            ! Calculate integrated intensity UIID

            IF (WIDE_BAND_MODEL .OR. WSGG_MODEL) THEN
               UIID(:,:,:,IBND) = UIID(:,:,:,IBND) + WEIGH_CYL*RSA(N)*IL
            ELSE
               UIID(:,:,:,ANGLE_INC_COUNTER) = UIID(:,:,:,ANGLE_INC_COUNTER) + WEIGH_CYL*RSA(N)*IL
            ENDIF

            ! Interpolate boundary intensities onto other meshes.
            ! IL_S is an array holding the intensities IL of cells just outside of mesh NOM.

            INTERPOLATION_LOOP: DO NOM=1,NMESHES
               IF (NM==NOM) CYCLE INTERPOLATION_LOOP
               M2=>OMESH(NOM)
               IF (M2%NIC_S==0) CYCLE INTERPOLATION_LOOP
               OTHER_WALL_LOOP: DO LL=1,M2%NIC_S
                  M2%IL_S(LL,N,IBND) = IL(M2%IIO_S(LL),M2%JJO_S(LL),M2%KKO_S(LL))
               ENDDO OTHER_WALL_LOOP
            ENDDO INTERPOLATION_LOOP

            ! Compute projected intensity on particles with a specified ORIENTATION

            IF (ORIENTED_PARTICLES) THEN
               PARTICLE_RADIATION_LOOP: DO IP=1,NLP
                  LP => LAGRANGIAN_PARTICLE(IP)
                  LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
                  IF (LPC%N_ORIENTATION==0) CYCLE PARTICLE_RADIATION_LOOP
                  BC => BOUNDARY_COORD(LP%BC_INDEX)
                  TEMP_ORIENTATION(1:3) = ORIENTATION_VECTOR(1:3,LP%ORIENTATION_INDEX)
                  IF (LP%INIT_INDEX > 0) THEN
                     IN => INITIALIZATION(LP%INIT_INDEX)
                     IF (ANY(IN%ORIENTATION_RAMP_INDEX > 0)) THEN
                        TEMP_ORIENTATION(1) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(1))
                        TEMP_ORIENTATION(2) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(2))
                        TEMP_ORIENTATION(3) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(3))
                        TEMP_ORIENTATION = TEMP_ORIENTATION / &
                                           (SQRT(TEMP_ORIENTATION(1)**2+TEMP_ORIENTATION(2)**2+TEMP_ORIENTATION(3)**2) &
                                           +TWO_EPSILON_EB)
                     ENDIF
                  ENDIF
                  COS_DLO = -DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N))
                  IF (COS_DLO > COS_HALF_VIEW_ANGLE(LP%ORIENTATION_INDEX)) THEN
                     DLO = -(TEMP_ORIENTATION(1)*DLX(N) + TEMP_ORIENTATION(2)*DLY(N) + TEMP_ORIENTATION(3)*DLZ(N))
                     BR => BOUNDARY_RADIA(LP%BR_INDEX)
                     IF (LPC%MASSLESS_TARGET) THEN
                        BR%BAND(IBND)%ILW(N) = DLO * IL(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_FACTOR(LP%ORIENTATION_INDEX)
                        IF (N==NEAREST_RADIATION_ANGLE(LP%ORIENTATION_INDEX)) &
                           BR%IL(IBND) = IL(BC%IIG,BC%JJG,BC%KKG)
                     ELSE
                        ! IL_UP does not account for the absorption of radiation within the cell occupied by the particle
                        BR%BAND(IBND)%ILW(N) = DLO * IL_UP(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_FACTOR(LP%ORIENTATION_INDEX)
                     ENDIF
                  ENDIF
               ENDDO PARTICLE_RADIATION_LOOP
            ENDIF

            ! Save the intensities for the radiation file (RADF)

            DO NN=1,N_RADF
               RF => RAD_FILE(NN)
               DO K=RF%K1,RF%K2,RF%K_STEP
                  DO J=RF%J1,RF%J2,RF%J_STEP
                     DO I=RF%I1,RF%I2,RF%I_STEP
                        RF%IL_SAVE(I,J,K,N) = IL(I,J,K)
                     ENDDO
                  ENDDO
               ENDDO
            ENDDO

         ENDDO ANGLE_LOOP

      ENDDO UPDATE_LOOP

      ! Compute incoming flux on walls and particles
      DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
         WC => WALL(IW)
         IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE
         B1 => BOUNDARY_PROP1(WC%B1_INDEX)
         SF => SURFACE(WC%SURF_INDEX)
         IF (SF%EXTERNAL_FLUX > TWO_EPSILON_EB) THEN
            IF (ABS(T_BEGIN) <= SPACING(B1%T_IGN)) THEN
               TSI = T
            ELSE
               TSI = T - B1%T_IGN
            ENDIF
            EFLUX = EVALUATE_RAMP(TSI,SF%RAMP(TIME_EFLUX)%INDEX,TAU=SF%RAMP(TIME_EFLUX)%TAU)*SF%EXTERNAL_FLUX
         ELSE
            EFLUX = 0._EB
         ENDIF
         B1%Q_RAD_IN  = B1%Q_RAD_IN + B1%EMISSIVITY*(INRAD_W(IW)+BBFA*EFLUX)
      ENDDO

      DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
         CFA => CFACE(ICF)
         IF (CFACE(ICF)%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE
         B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
         SF => SURFACE(CFA%SURF_INDEX)
         IF (SF%EXTERNAL_FLUX > TWO_EPSILON_EB) THEN
            IF (ABS(T_BEGIN) <= SPACING(B1%T_IGN)) THEN
               TSI = T
            ELSE
               TSI = T - B1%T_IGN
            ENDIF
            EFLUX = EVALUATE_RAMP(TSI,SF%RAMP(TIME_EFLUX)%INDEX,TAU=SF%RAMP(TIME_EFLUX)%TAU)*SF%EXTERNAL_FLUX
         ELSE
            EFLUX = 0._EB
         ENDIF
         B1%Q_RAD_IN  = B1%Q_RAD_IN + B1%EMISSIVITY*(INRAD_F(ICF)+BBFA*EFLUX)
      ENDDO

   ENDIF INTENSITY_UPDATE

   ! Save source term for the energy equation (QR = -DIV Q)

   IF (WIDE_BAND_MODEL .OR. WSGG_MODEL) THEN
      QR = QR + KAPPA_GAS*UIID(:,:,:,IBND)-KFST4_GAS
      IF (STORE_RADIATION_TERMS) THEN
         RADIATION_EMISSION   = RADIATION_EMISSION + KFST4_GAS
         RADIATION_ABSORPTION = RADIATION_ABSORPTION + KAPPA_GAS*UIID(:,:,:,IBND)
      ENDIF
      IF (N_LP_ARRAY_INDICES>0 .AND. PARTICLES_EXISTED) QR_W = QR_W + KAPPA_PART*UIID(:,:,:,IBND) - KFST4_PART
   ENDIF

ENDDO BAND_LOOP

! Sum up intensities and compute incoming flux at open boundaries

IF (UPDATE_INTENSITY) THEN

   UII = SUM(UIID, DIM = 4)

   DO IW=1,N_EXTERNAL_WALL_CELLS
      WC => WALL(IW)
      IF (WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) CYCLE
      B1 => BOUNDARY_PROP1(WC%B1_INDEX)
      BR => BOUNDARY_RADIA(WC%BR_INDEX)
      B1%Q_RAD_IN = 0._EB
      DO IBND=1,NUMBER_SPECTRAL_BANDS
         B1%Q_RAD_IN  = B1%Q_RAD_IN + SUM(BR%BAND(IBND)%ILW(1:NUMBER_RADIATION_ANGLES))
      ENDDO
   ENDDO

ENDIF

! Save source term for the energy equation (QR = -DIV Q) for the one-band (gray gas) case.
! QR for wide-band model is saved elsewhere.

IF (.NOT. (WIDE_BAND_MODEL .OR. WSGG_MODEL)) THEN
   QR = KAPPA_GAS*UII - KFST4_GAS
   IF (N_LP_ARRAY_INDICES>0) QR_W = QR_W + KAPPA_PART*UII - KFST4_PART
   IF (STORE_RADIATION_TERMS) THEN
      RADIATION_EMISSION   = KFST4_GAS
      RADIATION_ABSORPTION = KAPPA_GAS*UII
   ENDIF
ENDIF

! Calculate the incoming radiative flux onto the solid particles

IF (SOLID_PARTICLES .AND. UPDATE_INTENSITY) THEN
   PARTICLE_LOOP: DO IP=1,NLP
      LP => LAGRANGIAN_PARTICLE(IP)
      LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
      IF (LPC%SOLID_PARTICLE .OR. LPC%MASSLESS_TARGET) THEN
         SF => SURFACE(LPC%SURF_INDEX)
         B1 => BOUNDARY_PROP1(LP%B1_INDEX)
         IF (SF%EXTERNAL_FLUX > TWO_EPSILON_EB) THEN
            IF (ABS(T_BEGIN) <= SPACING(B1%T_IGN)) THEN
               TSI = T
            ELSE
               TSI = T - B1%T_IGN
            ENDIF
            EFLUX = EVALUATE_RAMP(TSI,SF%RAMP(TIME_EFLUX)%INDEX,TAU=SF%RAMP(TIME_EFLUX)%TAU)*SF%EXTERNAL_FLUX
         ELSE
            EFLUX = 0._EB
         ENDIF
         IF (LP%ORIENTATION_INDEX>0) THEN
            BR => BOUNDARY_RADIA(LP%BR_INDEX)
            B1%Q_RAD_IN = 0._EB
            DO IBND=1,NUMBER_SPECTRAL_BANDS
               B1%Q_RAD_IN = B1%Q_RAD_IN + B1%EMISSIVITY * (WEIGH_CYL*SUM(BR%BAND(IBND)%ILW(1:NUMBER_RADIATION_ANGLES)) + EFLUX)
            ENDDO
         ELSE
            BC => BOUNDARY_COORD(LP%BC_INDEX)
            B1%Q_RAD_IN = B1%EMISSIVITY*(0.25_EB*UII(BC%IIG,BC%JJG,BC%KKG) + EFLUX)
         ENDIF
      ENDIF
   ENDDO PARTICLE_LOOP
ENDIF

DEALLOCATE(IJK_SLICE)

! Write out intensities to the radiation file (RADF)

IF (N_RADF>0 .AND. T>=RADF_CLOCK(RADF_COUNTER(NM))) THEN

   DO N=1,N_RADF
      RF => RAD_FILE(N)
      OPEN(LU_RADF(N,NM),FILE=FN_RADF(N,NM),FORM='FORMATTED',STATUS='OLD',POSITION='APPEND')
      WRITE(LU_RADF(N,NM),'(/A)') 'TIME'
      WRITE(LU_RADF(N,NM),'(F8.2)') T
      WRITE(LU_RADF(N,NM),'(/A)') 'INTENSITIES'
      WRITE(FORMT,'(A,I4,A)') '(',NUMBER_RADIATION_ANGLES+2,'F12.2)'
      DO K=RF%K1,RF%K2,RF%K_STEP
         DO J=RF%J1,RF%J2,RF%J_STEP
            DO I=RF%I1,RF%I2,RF%I_STEP
               WRITE(LU_RADF(N,NM),FORMT) TMP(I,J,K),0._EB,(RF%IL_SAVE(I,J,K,NN),NN=1,NUMBER_RADIATION_ANGLES)
            ENDDO
         ENDDO
      ENDDO
      CLOSE(LU_RADF(N,NM))
   ENDDO

   RADF_COUNTER(NM) = RADF_COUNTER(NM) + 1

ENDIF

! UPDATE_ALL_ANGLES is a one-time only logical. If needed again, it will be reset again.

UPDATE_ALL_ANGLES = .FALSE.

END SUBROUTINE RADIATION_FVM


!> \brief Add user-specified HRRPUV to the heat release rate term, Q
!> \param MODE Indicator of whether volumetric heat release rate, HRRPUV, is to be added to or subtracted from Q
!> \details If MODE=0, add the user-specified volumetric HRRPUV to the total, Q (W/m3)
!> \details If MODE=1, add the user-specified volumetric HRRPUV to the total, Q (W/m3), and add radiative emission
!> \details If MODE=2, subtract the user-specified volumetric HRRPUV from the total, Q (W/m3)

SUBROUTINE ADD_VOLUMETRIC_HEAT_SOURCE(MODE)

INTEGER, INTENT(IN) :: MODE
REAL(EB) :: TIME_RAMP_FACTOR,RR
INTEGER :: N,I,J,K,ICC,JCC
TYPE(INITIALIZATION_TYPE), POINTER :: IN

DO N=1,N_INIT
   IN => INITIALIZATION(N)
   IF (IN%HRRPUV<=0._EB) CYCLE
   TIME_RAMP_FACTOR=1._EB
   IF (IN%RAMP_Q_INDEX>0) TIME_RAMP_FACTOR = EVALUATE_RAMP(T,IN%RAMP_Q_INDEX)
   DO K=0,KBP1
      DO J=0,JBP1
         DO I=0,IBP1
            SELECT CASE(IN%SHAPE)
               CASE DEFAULT
                  IF (XC(I) < IN%X1 .OR. XC(I) > IN%X2 .OR. &
                      YC(J) < IN%Y1 .OR. YC(J) > IN%Y2 .OR. &
                      ZC(K) < IN%Z1 .OR. ZC(K) > IN%Z2) CYCLE
               CASE('CONE')
                  RR = MAX(0._EB,IN%RADIUS*(1._EB-(ZC(K)-IN%Z0)/IN%HEIGHT))
                  IF ((XC(I)-IN%X0)**2+(YC(J)-IN%Y0)**2>RR**2) CYCLE
            END SELECT
            IF (MODE==1) THEN
               Q(I,J,K) = Q(I,J,K) + TIME_RAMP_FACTOR*IN%HRRPUV*IN%VOLUME_ADJUST(NM)
               KFST4_GAS(I,J,K) = KFST4_GAS(I,J,K) + IN%CHI_R*TIME_RAMP_FACTOR*IN%HRRPUV*IN%VOLUME_ADJUST(NM)
            ELSEIF (MODE==2) THEN
               Q(I,J,K) = Q(I,J,K) - TIME_RAMP_FACTOR*IN%HRRPUV*IN%VOLUME_ADJUST(NM)
            ELSE
               Q(I,J,K) = Q(I,J,K) + TIME_RAMP_FACTOR*IN%HRRPUV*IN%VOLUME_ADJUST(NM)
            ENDIF
            ! handle HRRPUV initialization in cutcells
            ! note: CC not yet connected to radiant fraction correction
            CC_IBM_IF: IF (CC_IBM) THEN
               IF (CCVAR(I,J,K,CC_CGSC) == CC_SOLID) EXIT CC_IBM_IF
               IF (CCVAR(I,J,K,CC_IDCC) > 0) THEN ! we have a cutcell
                  ICC=CCVAR(I,J,K,CC_IDCC)
                  DO JCC=1,CUT_CELL(ICC)%NCELL
                     CUT_CELL(ICC)%Q(JCC) = CUT_CELL(ICC)%Q(JCC) + TIME_RAMP_FACTOR*IN%HRRPUV
                  ENDDO
               ENDIF
            ENDIF CC_IBM_IF
         ENDDO
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE ADD_VOLUMETRIC_HEAT_SOURCE

END SUBROUTINE COMPUTE_RADIATION


REAL(EB) FUNCTION BLACKBODY_FRACTION(L1,L2,TEMP)

! Calculates the fraction of black body radiation between wavelengths L1 and L2 (micron) in Temperature TEMP

USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D_UNIFORM
REAL(EB),INTENT(IN) :: L1,L2,TEMP
REAL(EB) :: LT1,LT2,BBFLOW,BBFHIGH

LT1    =   L1 * TEMP/LTSTEP
CALL INTERPOLATE1D_UNIFORM(LBOUND(BBFRAC,1),BBFRAC,LT1,BBFLOW)

LT2    =   L2 * TEMP/LTSTEP
CALL INTERPOLATE1D_UNIFORM(LBOUND(BBFRAC,1),BBFRAC,LT2,BBFHIGH)

BLACKBODY_FRACTION = BBFHIGH - BBFLOW

END FUNCTION BLACKBODY_FRACTION


!> \brief Compute the radiative absorption coefficient, KAPPA
!> \param Z_IN Array of tracked species
!> \param TMP Temperature (K)
!> \param IBND Radiation band

REAL(EB) FUNCTION GET_KAPPA(Z_IN,TMP,IBND)

! Returns the radiative absorption

USE PHYSICAL_FUNCTIONS, ONLY : GET_MOLECULAR_WEIGHT
REAL(EB), INTENT(IN) :: Z_IN(1:N_TRACKED_SPECIES),TMP
REAL(EB) :: KAPPA_TEMP,INT_FAC,KAPPA_SUM,SCALED_Y_RADCAL_SPECIES,MWA
INTEGER, INTENT(IN) :: IBND
INTEGER :: LBND,UBND,N,TYY

KAPPA_SUM = 0._EB

TYY = MAX(0 , MIN(N_KAPPA_T,INT((TMP - RTMPMIN) * TYY_FAC)))

CALL GET_MOLECULAR_WEIGHT(Z_IN,MWA)

DO N = 1, N_RADCAL_ARRAY_SIZE
   SCALED_Y_RADCAL_SPECIES = DOT_PRODUCT(Z2RADCAL_SPECIES(N,:),Z_IN)
   IF (SCALED_Y_RADCAL_SPECIES<TWO_EPSILON_EB) CYCLE
   IF (RADCAL_SPECIES_INDEX(N)==16) THEN
      INT_FAC = MAX(0._EB,SCALED_Y_RADCAL_SPECIES)**0.25_EB
   ELSE
      INT_FAC = MAX(0._EB,SCALED_Y_RADCAL_SPECIES*MWA)**0.25_EB
   ENDIF
   LBND = INT(INT_FAC)
   INT_FAC = INT_FAC - LBND
   LBND = MIN(LBND,N_KAPPA_Y)
   UBND = MIN(LBND+1,N_KAPPA_Y)
   KAPPA_TEMP = RADCAL_SPECIES2KAPPA(N,LBND,TYY,IBND)
   KAPPA_SUM  = KAPPA_SUM + KAPPA_TEMP + INT_FAC*(RADCAL_SPECIES2KAPPA(N,UBND,TYY,IBND)-KAPPA_TEMP)
ENDDO

GET_KAPPA = KAPPA_SUM

END FUNCTION GET_KAPPA


!==================================================================================
!Function to compute the absorption coefficient according to Bordbar et al. (2014)
!==================================================================================
REAL(EB) FUNCTION KAPPA_WSGG(X_H2O, X_CO2, MOL_RATIO,PARTIAL_PRESSURE,JWSGG)
USE WSGG_ARRAYS

INTEGER, INTENT(IN) :: JWSGG
INTEGER :: NN
REAL(EB), INTENT(IN) :: X_H2O, X_CO2, MOL_RATIO,PARTIAL_PRESSURE
REAL(EB) :: SUM_KAPPA

! If no CO2 or H2O, return zero

IF ((X_H2O<=TWO_EPSILON_EB) .AND. (X_CO2<=TWO_EPSILON_EB)) THEN
   KAPPA_WSGG = 0._EB
   RETURN
ENDIF

!-------------------------------------------------------------------------
!Compute the absorption coefficient within three intervals of molar ratio
!-------------------------------------------------------------------------

IF (MOL_RATIO < 0.01_EB) THEN               !Only CO2
   SUM_KAPPA = WSGG_KAPPAP1_ARRAY(JWSGG)

ELSEIF (MOL_RATIO > 4._EB) THEN               !Only H2O
   SUM_KAPPA = WSGG_KAPPAP2_ARRAY(JWSGG)

ELSE                                          !CO2-H2O mixture
   IF (JWSGG==5) THEN
      SUM_KAPPA = 0._EB
   !Pressure-based absorption coefficient of the gray gases
   ELSE
      SUM_KAPPA = 0._EB
      DO NN = 0,4
         SUM_KAPPA = SUM_KAPPA + WSGG_D_ARRAY(JWSGG,NN)*(MOL_RATIO**REAL(NN,EB))
      ENDDO
   ENDIF
ENDIF

KAPPA_WSGG = SUM_KAPPA*PARTIAL_PRESSURE

END FUNCTION KAPPA_WSGG

!===================================================================================
!Function to compute the temperature coefficient according to Bordbar et al. (2014)
!===================================================================================
REAL(EB) RECURSIVE FUNCTION A_WSGG(TTMP,MOL_RATIO,JWSGG) &
   RESULT(A_FUNC_RES)
USE WSGG_ARRAYS

INTEGER,INTENT(IN) :: JWSGG
INTEGER :: MM,NN
REAL(EB),INTENT(IN) :: MOL_RATIO,TTMP
REAL(EB) :: TREF,SUM_A,SUM_B,SUM_C

!------------------------
!Parameters of the model
!------------------------
TREF = 1200._EB ! Reference temperature

!------------------------------------------------------------------
!Computing the temperature coefficient for the transparent windows
!------------------------------------------------------------------
IF (JWSGG==5) THEN
   SUM_A = 0._EB
   DO MM=1,4
      SUM_A = SUM_A + A_WSGG(TTMP,MOL_RATIO,MM)
   ENDDO
   A_FUNC_RES = 1._EB - SUM_A

ELSE
!---------------------------------------------------------
!Computing the temperature coefficient for the gray gases
!(within three intervals of molar ratio)
!---------------------------------------------------------
   IF (MOL_RATIO < 0.01_EB) THEN               !Only CO2
      !Computing the polynomial
      SUM_B = 0._EB
      DO MM=0,4
         SUM_B = SUM_B + WSGG_B1_ARRAY(JWSGG,MM)*(TTMP/TREF)**(REAL(MM,EB))
      ENDDO

   ELSEIF (MOL_RATIO > 4._EB) THEN               !Only H2O
      !Computing the polynomial
      SUM_B = 0._EB
      DO MM=0,4
         SUM_B = SUM_B + WSGG_B2_ARRAY(JWSGG,MM)*(TTMP/TREF)**(REAL(MM,EB))
      ENDDO

   ELSE
      !Computing the polynomials
      SUM_B = 0._EB
      DO MM=0,4
         SUM_C = 0._EB
         DO NN=0,4
            SUM_C = SUM_C + WSGG_C_ARRAY(JWSGG,MM,NN)*MOL_RATIO**(REAL(NN,EB))
         ENDDO
         SUM_B = SUM_B + SUM_C*(TTMP/TREF)**(REAL(MM,EB))
      ENDDO
   ENDIF
   A_FUNC_RES = SUM_B

ENDIF

END FUNCTION A_WSGG


!====================================================
!Function to compute the gray absorption coefficient
!of soot (same function as the one used by Fluent)
!====================================================
REAL(EB) FUNCTION KAPPA_SOOT(SOOT_MASS_CONCENTRATION,TTMP)

REAL(EB),INTENT(IN) :: SOOT_MASS_CONCENTRATION,TTMP

   !KAPPA_SOOT = 1232.4_EB*SOOT_MASS_CONCENTRATION*(1._EB+4.8E-4_EB*(TTMP-2000._EB))
   KAPPA_SOOT = SOOT_MASS_CONCENTRATION*(1232.4_EB+0.591552_EB*(TTMP-2000._EB))

END FUNCTION KAPPA_SOOT


END MODULE RAD




